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Executive  Summary  '  /  >.r 

/- 

Herein  we  address  the  question  of  whethe^the  techniques  of  nonlinear  dynamics 
system  theory -can  be  usefully  applied  in  a  physical  oceanography  context.  The  tentative 
answer  to  this  question  is*  yes,  the  dynamics  of  a  numerically  generated  one-dimensional 
sea  surface  is  shown  to  take  place  on  a  low-dimensional  attractor.  Further,  the  dimension 
of  the  attractor  is  fractional  (noninteger)  and  therefore  the  trajectory  describing  the  sur¬ 
face  evolution  is  chaotic. 

The  report  is  segmented  into  four  sections.  In  the  first  section  we  review  the  tradi¬ 
tional  wisdom  of  the  mathematical  modeling  of  waves  on  the  sea  surface.  In  the  second 
section  a  mini-review  of  nonlinear  dynamics  is  presented,  in  which  the  basic  concepts  of 
importance  in  understanding  the  influence  of  nonlinearities  on  the  evolution  of  a  system 
are  discussed  in  a  straightforward  way.  How  these  concepts  have  been  applied  in  geo¬ 
physical  context,  including  the  ocean  surface,  is  discussed  in  Section  3.  In  particular  it  is 
shown  that  both  climate  and  weather  have  chaotic  attractors  using  a  technique  that  allows 
one  to  reconstruct  the  attractor  directly  from  observational  data.  The  properties  of  the 
sea  surface  modeled  as  a  fractal  surface  are  also  discussed. 

In  Section  4  some  original  research  is  presented  in  which  the  attractor  reconstruc¬ 
tion  technique  is  applied  to  a  numerically  generated  one-dimensional  sea  surface  having 
a  Phillips  spectrum  of  waves,  In  these  preliminary  calculations  we  find  that,  like  the  cli¬ 
mate  and  weather  attractors ,^he  water  wave  attractor  has  a  low-order  fractional  dimen¬ 
sion.  This  implies  that  as  few  as  five  or  six  degrees  of  freedom  may  be  sufficient  to 
describe  the  dynamics  of  a  surface  that  required  512  degrees  of  freedom  to  numerically 
generate.  The  implications  of  this  for  such  problems  as  data  processing  are  still  being 
considered  and  one  discussed  inpart  in  Section  5.  We  emphasize  that  these  results  are 
preliminary  and  may  require  substantial  investigation  before  they  can  be  substantiated. 
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1.  Introduction 

It  is  well  known  that  the  geophysical  data  sets  seem  almost  capricious  in  the ir  varia¬ 
bility.  Indeed  it  is  this  mercurial  nature  of  geophysical  flows  that  makes  the  prediction  of 
weather  patterns,  or  sea  states,  from  deterministic  primitive  equations  so  uncertain.  One 
strategy  for  taking  this  changeability  into  account  has  been  to  introduce  noise  directly 
into  the  equations  of  motion,  thereby  replacing  the  deterministic  equations  of  motion  by 
stochastic  differential  equations  [Landau  and  Lifshitz,  1959].  From  these  stochastic 
equations  it  is  then  possible  to  construct  transport  equations  for  the  evolution  of  average 
flow  properties,  sucn  as  the  mean  velocity,  the  average  energy,  etc..  In  physical  oceanog¬ 
raphy,  at  least,  this  technique  has  not  been  very  successful  in  describing  the  evolution  of 
the  energy  spectral  density  of  gravity  waves  [see  eg..  West,  1981],  In  large  part  this 
failure  arises  because  such  transport  equations  do  not  properly  take  into  account  such 
strongly  nonlinear  effects  as  wave  breaking  and  long  wave/short  wave  interactions.  In 
this  report  we  wish  to  focus  on  certain  techniques  which  may  be  able  to  determine  the 
influence  of  strong  nonlinear  interactions  on  the  evolution  of  water  waves,  as  well  as  on 
weather  patterns  and  climate.  These  techniques  have  been  developed  in  the  area  of  non¬ 
linear  systems  theory. 

In  Section  2  of  this  report  we  present  a  mini-review  of  nonlinear  dynamics  which 
may  be  skipped  over  by  the  experts,  but  which  the  novice  might  find  useful.  In  this  sec¬ 
tion  one  will  discover  that  dynamic  systems  theory  emerged  from  a  fusion  of  two  classi¬ 
cal  areas  of  mathematics,  topology  and  the  theory  of  differential  equations.  Its  impor¬ 
tance  to  the  experimental  and  observational  sciences  lies  in  its  capacity  to  quantitatively 
characterize  complex  dynamical  behavior.  Herein  we  review  how  dynamical  systems 
theory  is  applied  in  various  geophysical  contexts.  One  way  in  which  it  is  applied  is  in  the 
construction  of  simple  dynamical  models  that  give  rise  to  solutions  that  resemble  the 
observed  time  series  data.  Another  way  in  which  it  is  applied  is  through  the  development 
of  data  processing  algorithms  that  capture  the  essential  features  of  the  dynamics  of  the 
system,  such  as  its  degree  of  irregularity  and  the  structure  of  the  attractor*  on  which  the 
system’s  dynamics  takes  place  [cf.  Section  2],  It  is  obvious  that  the  theory  of  differential 
equations  is  useful  because  it  enables  one  to  construct  the  dynamic  equations  that 


All  technical  terms,  from  nonlinear  dynamics,  used  in  the  Introduction  are  defined  in  Section  2. 
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describe  the  evolution  of  the  geophysical  system  of  interest.  Topology  is  of  value  here 
because  it  allows  us  to  determine  the  unique  geometrical  properties  of  the  resulting 
dynamic  attractors.  The  degree  of  irregularity  or  randomness  of  calculated  time  series  is 
closely  related  to  the  geometrical  structure  of  the  underlying  attractor. 

Mandelbrot  (1977,  1982),  the  father  of  fractals,  argues  that  we  need  a  new  kind  of 
geometry  to  study  such  structures.  Euclidean  geometry  is  concerned  with  the  under¬ 
standing  of  straight  lines  and  regular  forms,  and  it  is  usually  assumed,  in  geophysical  and 
oceanography,  that  the  world  consists  of  continuous,  smooth  curves  in  spaces  of  integer 
dimension.  When  we  look  at  billowing  clouds,  trees  of  all  kinds,  breaking  waves  and 
coastlines,  we  observe  that  the  notions  of  classical  geometry  are  inadequate  to  describe 
them.  Detail  does  not  become  less  important  as  regions  of  these  various  structures  are 
magnified,  but  perversely  more  and  more  detail  is  revealed  at  each  level  of  magnification. 
The  rich  texture  of  these  structures  is  characteristic  of  fractals  [Mandelbrot,  1977,  1982]. 
In  Section  3  we  show  that  a  fractal  structure  is  not  smooth  and  homogeneous,  and  that  the 
smaller-scale  structure  is  similar  to  the  large-scale  form.  The  aspect  of  such  structures 
that  makes  them  different  from  what  we  usually  experience  is  that  there  is  no  characteris¬ 
tic  length  scale. 

But  it  is  not  only  static  structures  that  have  fractal  properties  but  dynamical 
processes  as  well.  In  Section  3  we  also  review  how  the  fractal  concept  has  been  applied 
to  time  series  in  some  geophysical  applications  and  in  Section  4  some  preliminary  appli¬ 
cations  are  made  to  a  spectrum  of  deep  water  gravity  waves.  Because  this  latter  applica¬ 
tion  is  particularly  important  to  us,  let  us  now  review  some  of  the  traditional  wisdom 
regarding  such  wave  fields. 

The  generation,  growth,  propagation  and  eventual  dissipation  of  waves  on  the  ocean 
surface  have,  over  the  past  quarter  century,  been  described  in  terms  of  weakly  interact¬ 
ing,  nonlinear  waves  [see  eg.,  Phillips,  1966;  West,  1981].  The  equations  of  motion  for 
the  nonlinear  wave  field  are  determined  by  a  Hamiltonian  [Zakhorov,  1968;  Miles,  1977] 
for  an  incompressible,  inviscid,  irrotational  fluid.  To  describe  the  process  of  wave  gen¬ 
eration,  evolution,  and  the  subsequent  development  of  wave  instabilities,  it  is  often  con¬ 
venient  to  express  the  observables  at  the  ocean  surface  as  series  expansions  in  the  eigen¬ 
function  of  the  linearized  system.  For  water  waves  these  eigenfunctions  are  simple  sines 
and  cosines  and  the  eigenfunction  expansion  is  merely  a  Fourier  decomposition  of  the 
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surface.  The  expansion  coefficients  are  interpreted  as  the  constant  amplitudes  of 
independent  waves  in  a  linear  wave  field,  and  as  the  time  varying  amplitudes  in  a  non¬ 
linear  wave  field.  The  Hamiltonian  for  this  system  can  be  expressed  as  a  series  in  which 
the  nonlinear  terms  appear  as  products  of  the  mode  amplitudes,  see  eg.  West  (1981). 
These  nonlinear  interactions,  interpreted  as  the  scattering  or  coupling  of  the  once  linear 
waves,  induce  variations  in  both  the  amplitudes  and  phases  cf  the  linear  modes  in  the 
equations  of  motion. 

The  sea  surface  gravity  wave  field  is  described  as  a  conservative  Hamiltonian  sys¬ 
tem  and  Hamilton’s  equations  of  motion  provides  a  deterministic  description  of  its  evolu¬ 
tion  [Hasselmann,  1968;  Broer,  1974;  Watson  and  West,  1975;  Miles,  1977;  West,  1981], 
If  we  assume  that  this  field  is  well  represented  by  N  degrees  of  freedom,  where  N  may 
be  large  but  finite  (a  field  has  an  infinite  number  of  degrees  of  freedom),  then  the  system 
can  be  represented  by  2 N  coupled,  deterministic,  nonlinear  rate  equations  for  the  mode 
amplitudes.  Moser  (1973)  separates  the  interactions  into  resonant  and  nonresonant 
groups.  The  nonresonant  interactions  provide  for  a  stable  evolution  of  the  wave  field, 
whereas  the  resonant  interactions  produce  instabilities.  The  existence  of  resonances  in 
water  wave  dynamics  was  explicitly  pointed  out  by  Phillips  (1960).  He  demonstrated 
that  just  as  for  resonances  in  linear  systems,  the  nonlinear  resonances  in  wave/wave 
interactions  produce  an  initial  secular  growth  of  new  waves.  Benny  (1962)  extended 
these  arguments  to  show  how  the  nonlinear  interactions  also  saturate,  thereby  quenching 
the  apparent  instability.  Chirikov  (1979)  pointed  out  that  the  instabilities  generated  by 
such  nonlinear  resonances  are  always  bounded,  unlike  linear  resonances  which  grow 
without  limit.  This  bounding  of  the  instability  is  produced  by  the  nonlinear  dispersion 
relation,  i.e.,  the  amplitude  dependence  of  the  frequencies  in  the  equations  of  motion. 
Weak  nonlinearities,  therefore,  act  to  stabilize  the  water  wave  system  and  inhibit  explo¬ 
sive  instabilities. 

The  concept  of  dominant  resonant  interactions  have  formed  the  basis  for  many 
analytic-numeric  calculations  of  the  properties  of  surface  water  waves  [Hasselmann, 
1962,  1963a,  b;  Zakhorov,  1968;  Yuen  and  Lake,  1982;  Peregrine,  1983;  Bryant,  1984], 
The  evolution  of  capillary  and  gravity-capillary  waves  has  been  described  by  mode  rate 
equations  with  cubic  nonlinearities,  i.e.,  four-wave  interactions.  The  nonresonant 
interactions  had  been  thought  to  be  unimportant  in  describing  the  evolution  of  the  wave 
field;  however,  Watson  and  West  (1975),  using  perturbation  theory,  and  West  (1981), 
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using  canonical  transformations  of  variables,  demonstrated  that  the  nonresonant  terms 
contribute  substantially  to  the  strength  of  the  resonant  interactions  among  the  surface 
waves  and  are  therefore  not  negligible.  Thus,  in  the  numerical  results  discussed  in  Sec¬ 
tion  4,  we  include  both  resonant  and  nonresonant  effects. 

The  motion  of  the  free  surface  of  an  incompressible,  irrotational,  inviscid  fluid 
£  ( x,  t )  is  determined  by  Bernoulli’s  equation 

-jUt>  +  %v-v  +  g C  =  0  ,z  -  C  ,  (1-1) 

with  the  fluid  velocity  v  ( x,  z ,  t )  given  by  the  gradient  of  flic  velocity  potential  <J>  ( x,  z ,  t ), 

v  =  V<j>  (1.2) 

and  the  vertical  motion  of  the  free  surface  by 

■—  £(x,i)  +  v-V£(x,r)  =  =  W  ,z  =S  .  (1.3) 


In  the  interior  of  the  incompressible  fluid, 

Vv  =  0 

so  that  the  velocity  potential  satisfies  Laplace’s  equation 

V2<}>  =  0  . 

These  equations  can  be  converted  to  equations  at  the  free  surface,  where 

<\>s  (x,r)  s<}>[x,z  =  C(x,0.*l 
with  the  result  [Watson  and  West,  1975], 


+//z(VJ<DJ)2  +  g^  =  '/z[l  +  (V,C)2]W2  ,  (1.7) 

+  ,  (1,8) 

where  is  the  horizontal  gradient  operator.  The  numerical  integration  of  (1.7)  and 
(1.8)  is  based  on  the  method  of  Watson  and  West  (1975)  and  is  done  by  taking  the  indi¬ 
cated  products  of  the  field  quantities  in  configuration  space;  fast  Fourier  transforming 
(FFT)  the  configuration  equations  md  time  incrementing  the  transformed  equations  to 
obtain  the  components  of  the  appropriate  field  variables,  then  transforming  (FFT)  back  to 
configuration  space  to  again  evaluate  the  nonlinear  products  and  start  the  process  again. 
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This  numerical  procedure,  developed  by  West,  Brueckner,  Janda,  Milder  and  Milton 
(1987),  does  not  distinguish  between  resonant  and  nonresonant  interactions,  so  that  it 
includes  terms  explicitly  neglected  or  approximated  by  most  calculations  involving  the 
use  of  nonlinear  mode  rate  equations. 

The  question  of  interest  here  is  whether  the  sea  surface  described  by  these  equa¬ 
tions,  and  therefore  explicitly  determined  by  the  evolution  of  2 N  degrees  of  freedom  can 
also  be  described  by  a  low-dimensional  dynamic  attractor.  This  question  is  addressed  in 
Section  4,  where  the  techniques  discussed  in  Section  2  and  3  are  applied  to  the  time 
series  generated  by  the  numerical  code.  The  tentative  answer  to  the  question  is,  yes,  the 
dynamics  of  the  one-dimensional  sea  surface  does  take  place  on  a  low-dimensional 
attractor.  Further,  the  attractor  seems  to  have  a  fractional  (fractal)  dimension  and  there¬ 
fore  the  trajectory  describing  the  evolution  is  chaotic.  We  emphasize  that  these  results 
are  preliminary  and  may  require  substantial  investigation  before  they  can  be  substan¬ 
tiated. 
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2.  Mini-Review  of  Nonlinear  Dynamics 

In  this  section  we  sketch  the  concepts  from  nonlinear  dynamics  that  will  be  the  most 
useful  in  attempting  to  determine  if  the  surface,  or  interior,  of  the  ocean  can  be  usefully 
modeled  as  a  low-dimensional  attractor.  One  of  the  most  fruitful  and  brilliant  ideas  of  the 
second  half  of  the  1600’s  was  the  idea  that  the  concept  of  a  function  and  the  geometric 
representation  of  a  curve  are  related.  Geometrically  the  notion  of  a  linear  relation 
between  two  quantities  implies  that  if  a  graph  is  constructed  with  the  ordinate  denoting 
the  values  of  one  variable  and  the  abscissa  denoting  the  values  of  the  other  variable,  then 
the  relation  in  question  appears  as  a  straight  line.  In  general,  the  graph  or  curve  lies  in 
the  “space”  defined  by  the  independent  coordinate  axes.  In  a  dynamic  system  the  coor¬ 
dinate  axis  are  defined  by  the  possible  values  of  the  independent  dynamic  variables. 

The  state  of  a  given  system  is  defined  by  a  point  in  the  above  space,  often  called 
either  the  state  space  or  phase  space  of  the  system.  As  time  moves  on  the  point  traces  out 
a  curve,  called  an  orbit  or  trajectory,  which  describes  the  history  of  the  system’s  evolu¬ 
tion.  This  geometrical  representation  of  dynamics  is  one  of  the  most  useful  tools  in 
dynamic  systems  theory  for  analyzing  the  time-dependent  properties  of  nonlinear  sys¬ 
tems.  By  nonlinear  we  mean  the  output  of  a  system  is  not  proportional  to  the  input.  One 
implication  of  this  is  the  following:  If  the  system  is  linear,  than  two  trajectories  initiated 
at  nearby  points  in  phase  space  evolve  in  close  proximity,  so  that  at  any  point  in  future 
time  the  two  trajectories  (and  therefore  the  states  of  the  system  they  represent)  are  also 
near  to  one  another.  If  the  system  is  nonlinear  then  two  such  trajectories  could  diverge 
from  one  another  and  at  subsequent  times  the  two  trajectories  could  be  arbitrarily  far 
apart,  i.e.,  the  distance  between  the  orbits  does  not  evolve  in  a  proportionate  way.  Of 
course  this  need  not  necessarily  happen  in  a  nonlinear  system,  it  is  a  question  of  stability. 

This  brings  us  back  to  our  recurrent  example  of  the  weather  and  to  the  question  of 
its  predictability.  Its  elusive  nature  has  only  recently  come  into  sharper  focus  and  made 
clear  two  distinct  views  of  the  character  of  the  evolution  of  deterministic  dynamic  sys¬ 
tems.  These  views  were  articulated  by  their  respective  proponents,  Laplace  and 
Poincar^,  writing  nearly  100  years  apart. 
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Laplace  (1825)  believed  in  strict  determinism  and  to  his  mind  this  implied  complete 
predictability.  Uncertainty  for  him  is  a  consequence  of  imprecise  knowledge,  so  that 
probability  theory  is  necessitated  by  incomplete  and  imperfect  observations.  This  view  is 
clearly  expressed  in  the  quotation: 


‘  ‘The  present  state  of  the  system  of  nature  is  evidently  a  consequence  of  what  it  was 
in  the  preceding  moment,  and  if  we  conceive  of  an  intelligence  which  at  a  given  instant 
comprehends  all  the  relations  of  the  entities  of  this  universe,  it  could  state  the  respective 
positions,  motions,  and  general  affects  of  all  these  entities  at  any  time  in  the  past  or 
future.” 

“Physical  astronomy,  the  branch  of  knowledge  which  does  the  greatest  honor  to  the 
human  mind,  gives  us  an  idea,  albeit  imperfect,  of  what  such  an  intelligence  would  be. 
The  simplicity  of  the  law  by  which  the  celestial  bodies  move,  and  the  relations  of  their 
masses  and  distances,  permit  analysis  to  follow  their  motions  up  to  a  certain  point;  and  in 
order  to  determine  the  state  of  the  system  of  these  great  bodies  in  past  or  future  centuries, 
it  suffices  for  the  mathematician  that  their  position  and  their  velocity  be  given  by  obser¬ 
vation  for  any  moment  in  time.  Man  owes  that  advantage  to  the  power  of  the  instrument 
he  employs,  and  to  the  small  number  of  relations  that  it  embraces  in  its  calculations.  But 
ignorance  of  the  different  causes  involved  in  the  productions  of  events,  as  well  as  their 
complexity,  taken  together  with  the  imperfection  of  analysis,  prevents  our  reaching  the 
same  certainty  about  the  vast  majority  of  phenomena.  Thus  there  are  things  that  are  unc¬ 
ertain  for  us,  things  more  or  less  probable,  and  we  seek  to  compensate  for  the  impossibil¬ 
ity  of  knowing  them  by  determining  their  different  degrees  of  likelihood.  So  it  is  that  we 
owe  to  the  weakness  of  the  human  mind  one  of  the  most  delicate  and  ingenious  of 
mathematical  theories,  the  science  of  chance  or  probability.” 


Poincard  1906  on  the  other  hand  sees  an  intrinsic  inability  to  make  predictions  due 
to  a  sensitive  dependence  of  the  evolution  of  the  system  on  the  initial  state  of  the  system. 
He  expressed  this  in  the  following  way: 

“A  very  small  cause  which  escapes  our  notice  determines  a  considerable  effect  that 
we  cannot  fail  to  see,  and  then  we  say  that  the  effect  is  due  to  chance.  If  we  knew 
exactly  the  laws  of  nature  and  the  situation  of  the  universe  at  the  initial  moment,  we 
could  predict  exactly  the  situation  of  that  same  universe  at  a  succeeding  moment.  But 
even  if  it  were  the  case  that  the  natural  laws  had  no  longer  any  secret  for  us,  we  could 
still  only  know  the  initial  situation  approximately.  If  that  enabled  us  to  predict  the 
succeeding  situation  with  the  same  approximation,  that  is  all  we  require,  and  we  should 
say  that  the  phenomenon  had  been  predicted,  that  it  is  governed  by  laws.  But  it  is  not 
always  so;  it  may  happen  that  small  differences  in  the  initial  conditions  produce  very 
great  ones  in  the  final  phenomena.  A  small  error  in  the  former  will  produce  an  enormous 
error  in  the  latter.  Prediction  becomes  impossible,  and  we  have  the  fortuitous 
phenomenon.” 
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Thus  we  can  clearly  distinguish  between  the  “traditional”  view  of  Laplace  and  the 
“modem”  view  of  Poincard  The  latter  is  considered  modem  because  we  have  now 
discovered  that  deterministic  systems  with  only  a  few  degrees  of  freedom  can  generate 
aperiodic  behavior  that  for  many  purposes  is  indistinguishable  from  random  fluctuations. 
We  emphasize  that  the  random  aspect  is  fundamental  to  the  system  dynamics  and  gather¬ 
ing  more  information  will  not  reduce  the  degree  of  uncertainty.  Randomness  generated 
in  this  way  is  now  called  chaos  [see  eg.  Crutchfield  et  al.,  1987].  To  understand  chaos 
we  need  to  discuss  the  dynamics  of  nonlinear  systems. 

We  introduced  the  notion  of  a  phase  space  and  a  trajectory  to  describe  the  dynamics 
of  a  system.  Each  choice  of  an  initial  state  for  the  system  produces  a  different  trajectory. 
If,  however,  there  is  a  limiting  set  in  phase  space  to  which  all  trajectories  are  drawn  after 
a  sufficiently  long  time,  we  say  that  the  system  dynamics  are  asymptotically  described  by 
an  attractor.  The  attractor  is  the  geometric  limiting  set  on  which  all  the  trajectories  even¬ 
tually  find  themselves,  i.e.,  the  set  of  points  in  phase  space  to  which  the  trajectories  are 
attracted.  Attractors  come  in  many  shapes  and  sizes,  but  they  all  have  the  property  of 
occupying  a  finite  volume  of  phase  space.  As  a  system  evolves  it  sweeps  through  the 
attractor,  going  through  some  regions  rather  rapidly  and  others  quite  slowly,  but  always 
staying  on  the  attractor.  Whether  or  not  the  system  is  chaotic  is  determined  by  how  two 
initially  adjacent  trajectories  cover  the  attractor  over  time.  As  Poincar£  stated,  if  a  small 
change  in  the  initial  separation  of  trajectories  (error)  produces  an  enormous  change  in 
their  final  separation  (error),  then  the  evolution  is  unpredictable.  One  question  is  how 
this  growing  separation,  indicative  of  chaos,  is  accomplished  on  an  attractor  of  finite  size. 
The  answer  has  to  do  with  the  layered  structure  of  the  attractor  necessary  for  it  to  be 
chaotic. 

Rossler  (1976)  described  chaos  as  resulting  from  the  geometric  operations  of 
stretching  and  folding.  Two  initially  nearby  orbits  cannot  rapidly  separate  forever  on  an 
attractor  of  finite  size,  therefore  the  attractor  must  eventually  fold  over  onto  itself.  Once 
folded  the  attractor  is  again  stretched  and  folded  again.  This  process  is  repeated  over  and 
over  yielding  an  attractor  structure  with  an  infinite  number  of  layers  to  be  traversed  by 
the  various  trajectories.  The  infinite  richness  of  the  attractor  structure  affords  ample 
opportunity  for  trajectories  to  diverge  and  follow  increasingly  different  paths.  The  finite 
size  of  the  attractor  insures  that  these  diverging  trajectories  will  eventually  pass  close  to 
one  another  again,  albeit  on  different  layers  of  the  attractor.  Crutchfield  et  al.  (1987) 
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visualize  thes  orbits  on  a  chaotic  attractor  as  being  shuffled  by  this  stretching  and  folding 
process,  much  as  a  deck  of  cards  is  shuffled  by  a  dealer.  Thus  the  randomness  of  the 
chaotic  orbits  is  a  consequence  of  this  shuffling  process.  This  operation  of  stretching  and 
folding  creates  folds  within  folds  ad  infinitum,  resulting  in  the  attractor  often  having  a 
fractal  structure  in  phase  space. 

The  degree  of  irregularity  (chaos)  of  the  dynamic  observable  is  closely  related  to  the 
geometrical  structure  of  the  underlying  attractor.  There  are  a  number  of  measures  of  the 
degree  of  chaos  of  these  attractors.  One  is  its  dimension-,  integer  values  of  the  dimension 
indicate  a  simple  attractor,  non-integer  dimension  indicates  a  chaotic  attractor  in  phase 
space. 

A  second  measure  of  the  degree  of  irregularity  generated  by  a  chaotic  attractor  is 
the  “entropy”  of  the  motion.  The  entropy  is  interpreted  by  Crutchfield  et  al.  (1987)  as 
the  average  rate  of  stretching  and  folding  of  the  attractor,  or  alternatively,  as  the  average 
rate  at  which  information  is  generated.  The  application  of  the  information  concept  in  the 
dynamic  systems  context  has  been  championed  by  Shaw  (1981,1984)  and  Nicolis 
(1985,1986).  One  can  view  the  preparation  of  the  initial  state  of  the  system  as  initializ¬ 
ing  a  certain  amount  of  information.  The  more  precisely  the  initial  state  can  be  specified, 
the  more  information  one  has.  This  corresponds  to  localizing  the  initial  state  of  the  sys¬ 
tem  in  phase  space,  the  amount  of  information  is  inversely  proportional  to  the  volume  of 
state  space  localized  by  measurement.  In  a  regular  attractor,  trajectories  initiated  in  a 
given  local  volume  stay  near  to  one  another  as  the  system  evolves,  so  the  initial  informa¬ 
tion  is  preserved  in  time  and  no  new  information  is  generated.  Thus  the  initial  informa¬ 
tion  can  be  used  to  predict  the  final  state  of  the  system.  In  a  chaotic  attractor  the  stretch¬ 
ing  and  folding  operations  smear  out  the  initial  volume,  thereby  destroying  the  initial 
information  as  the  system  evolves  and  the  dynamics  create  new  information.  Thus  the 
initial  uncertainty  in  the  specification  of  the  system  is  eventually  smeared  out  over  the 
entire  attractor  and  all  predictive  power  is  lost,  ie.,  all  causal  connection  between  the 
present  and  the  future  is  lost. 

Let  us  denote  the  region  of  phase  space  as  initially  occupied  by  Vi  (initial  volume) 
and  the  final  region  by  Vy.  The  change  in  the  observable  information  I  is  then  [Shaw, 
1981;  Nicolis  and  Tsuda,  1985] 
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The  rate  of  information  creation  or  dissipation  is  given  by 

dJ_  _  1  dV 
dt  V  dt 


(2.0.2) 


In  non-chaotic  systems,  the  sensitivity  of  the  flow  to  the  initial  conditions  grows  with 
time  at  most  as  a  polynomial,  eg.,  let  co(t )  be.the  number  of  distinguishable  states  so  that 


co(r)  -  tn 


since  VJVt  =  (Oy/co,-  we  have  [Shaw  19811 


dl_ _ n_ 

dt  t 


(2.0.3) 


(2.0.4) 


Thus  the  rate  of  information  generation  converges  to  zero  as  t  — »«>  and  the  final  state  is 
predictable  from  the  initial  information.  On  the  other  hand,  in  chaotic  systems  the  sensi¬ 
tivity  of  the  flow  to  the  initial  conditions  grow  exponentially  with  time, 


co(/)  -  e' 


so  that 


(2.0.5) 


(2.0.6) 


This  latter  system  is  therefore  a  continuous  source  of  information,  the  attractor  itself  gen¬ 
erates  the  information  independently  of  the  initial  conditions. 

The  final  measure  of  the  degree  of  chaos  associated  with  an  attractor  with  which  we 
will  be  concerned  is  the  set  of  Lyapunov  exponents.  These  exponents  quantify  the  aver¬ 
age  exponential  convergence  or  divergence  of  nearby  trajectories  in  the  phase  space  of 
the  dynamical  systems.  Wolf,  Swift,  Swinney  and  Vastano  (1985),  among  others, 
believe  the  spectrum  of  Lyapunov  exponents  provides  the  most  complete  qualitative  and 
quantitative  characterization  of  chaotic  behavior.  A  system  with  one  or  more  positive 
Lyapunov  exponents  is  defined  to  be  chaotic.  The  local  stability  properties  of  a  system 
are  determined  by  its  response  to  perturbations;  along  certain  directions  the  response  can 
be  stable  whereas  along  others  it  can  be  unstable.  If  we  consider  a  d-dimensional  sphere 
of  initial  conditions  and  follow  the  evolution  of  this  sphere  in  time,  then  in  some  direc¬ 
tions  the  sphere  will  contract,  whereas  in  others  it  will  expand,  thereby  forming  a  d- 
dimensional  ellipsoid.  Thus,  a  d-dimensional  system  can  be  characterized  by  d 
exponents  where  the  jth  Lyapunov  exponent  quantifies  the  expansion  or  contraction  of 
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the  flow  along  the  jth  ellipsoidal  principal  axis.  The  sum  of  the  Lyapunov  exponents  is 
the  average  divergence  of  the  flow,  which  for  a  dissipative  system  (possessing  an  attrac¬ 
tor)  must  always  be  negative. 

Consider  a  three  dimensional  phase  space  in  which  the  limiting  set  (the  attractor) 
can  be  characterized  by  the  triple  of  Lyapunov  exponents  (Xj,^,^).  The  qualitative 
behavior  of  the  attractor  can  be  specified  by  determining  the  signs  of  the  Lyapunov 
exponents  only,  ie.,  (signX^signy^jignXj).  As  shown  in  Figure  (2.0.1a)  the  triple  (-,-,-) 
corresponds  to  an  attracting  fixed  point.  In  each  of  the  three  directions  there  is  an 
exponential  contraction  of  trajectories,  so  that  no  matter  what  the  initial  state  of  the  sys¬ 
tem  it  will  eventually  wind  up  at  the  fixed  point.  This  fixed  point  need  not  be  the  origin, 
as  it  would  be  for  a  dissipative  linear  system,  but  can  be  anywhere  in  phase  space.  The 
arrows  shown  in  the  figure  do  not  necessarily  represent  trajectories  since  the  fixed  point 
can  be  approached  at  any  angle  by  the  evolving  nonlinear  system. 

An  attracting  limit  cycle  is  denoted  by  (0,-,-)  in  which  there  are  two  contracting 
directions  and  one  that  is  neutrally  stable.  In  Figure  (2.0.1b)  we  see  that  this  attractor 
resembles  the  orbit  of  a  harmonic  oscillator  with  a  particular  energy,  but  that  is  not  the 
case.  The  orbit  of  a  harmonic  oscillator  does  not  attract  points  from  off  the  orbit  onto 
itself.  On  the  other  hand  in  a  nonlinear  dynamical  system  an  orbit  has  a  basin  of  attrac¬ 
tion  so  that  all  systems  whose  initial  state  lies  in  the  basin  eventually  wind  up  on  the  limit 
cycle. 

The  triple  (0,0,-)  has  two  neutral  directions  and  one  that  is  contracting  so  that  the 
attractor  is  the  2-torus  depicted  in  Figure  (2.0.1c).  An  example  of  such  a  system  would 
be  two  coup  led  harmonic  oscillators,  where  the  two  positions  and  two  velocities  would 
describe  the  dynamics.  The  constant  energy  (no  dissipation)  reduces  the  number  of  vari¬ 
ables  in  this  coupled  system  to  three  so  that  the  system  is  described  by  the  two  constant 
radii  and  the  two  angles  locating  the  trajectory  on  the  surface  of  the  torus. 

Finally  (+,0,-)  corresponds  to  a  chaotic  attractor  in  which  the  trajectories  expand  in 
one  direction,  are  neutrally  stable  in  another  and  contracting  in  a  third.  In  order  for  the 
trajectories  to  continuously  expand  in  one  direction  and  yet  remain  on  a  finite  attractor, 

*  By  "flow''  we  mean  the  behavior  on  phase  space  of  a  bundle  of  trajectories  having  a  distribution  of  initial  conditions. 
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the  attractor  must  undergo  the  stretching  and  folding  operations  in  this  direction  as  dis¬ 
cussed  by  Rossler  (1978).  Figure  (2.0.  Id)  is  the  two-dimensional  projection  of  the 
Rossler  attractor. 

The  resolution  of  the  apparent  conflict  between  the  traditional  and  the  modem  view 
of  dynamic  systems  theory  as  presented  in  classical  mechanics  is  that  chaos  is  not  incon¬ 
sistent  with  the  traditional  notion  of  solving  deterministic  equations  of  evolution.  As 
Ford  (1987)  states: 

“...Determinism  means  that  Newtonian  orbits  exist  and  are  unique,  but  since 
existence- uniqueness  theorems  are  generally  nonconstructive,  they  assert  nothing  about 
the  character  of  the  Newtonian  orbits  they  define.  Specifically,  they  do  not  preclude  a 
Newtonian  orbit  from  passing  every  computable  test  for  randomness  of  being  humanly 
indistinguishable  from  a  realization  of  a  truly  random  process.  Thus,  popular  opinion  to 
the  contrary  notwithstanding,  there  is  absolutely  no  contradiction  in  the  term  “determin¬ 
istically  random.”  Indeed,  it  is  quite  reasonable  to  suggest  that  the  most  general 
definition  of  chaos  should  read:  chaos  means  deterministically  random...” 

2.1  Nonlinear  Oscillator 

In  the  physical  sciences  the  dynamics  of  a  system  are  determined  by  the  equations 
describing  how  the  physical  observables  change  in  time.  These  equations  are  obtained 
by  means  of  some  general  principle,  such  as  the  conservation  of  energy  and/or  momen¬ 
tum,  applied  to  the  system  of  interest.  The  appropriate  conservation  law  follows  from  a 
symmetry  of  the  system  which  determines  a  rule  by  which  the  system  evolves.  If  a  set  of 
circumstances  is  specified  by  an  N -component  vector  X  =  (X  j,X2, ..,  XN )  then  in  order  to 
predict  the  future  state  of  the  system  from  its  present  configuration,  we  must  specify  a 
rule  for  the  systems’  evolution.  In  the  physical  sciences  the  traditional  strategy  is  to  con¬ 
struct  a  set  of  differential  equations.  These  equations  are  obtained  by  considering  each 
component  of  the  system  to  be  a  function  of  time,  then  as  time  changes  so  too  do  the  cir¬ 
cumstances.  If  in  a  short  time  interval  A  t  we  can  associate  an  attendant  set  of  changes 
AX  =  (AXj, ...,  &XN)  as  determined  by  A X  =  F  (X,r)Ar  then  in  the  limit  Ar  — >0  one 
would  write  the  “equations  of  motion” 

-f  X(f)  =  F(X,0  (2.1.1) 

dt 

which  is  a  statement  about  the  evolution  of  the  system  in  time.  If  at  time  t  =  0  we 
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specify  the  components  X(0),  i.e.,  the  set  of  circumstances  characterizing  the  system,  and 
if  F(X,  t)  is  an  analytic  function  of  its  arguments,  then  the  evolution  of  the  system  is 
determined  by  direct  integration  of  the  equations  of  motion  away  from  the  initial  state. 

The  mathematicians  have  categorized  the  solutions  to  such  equations  for  the  sim¬ 
plest  kinds  of  systems.  One  way  to  describe  such  systems  is  by  means  of  geometric  con¬ 
structions  in  which  the  solution  to  an  equation  of  the  above  form  is  depicted  by  a  curve  in 
an  appropriate  space.  The  coordinate  axes  necessary  for  such  a  construction  are  the  con¬ 
tinuum  of  values  that  the  vector  X(/ )  can  assume,  each  axis  is  associated  with  one  com¬ 
ponent  of  the  vector  X;  this  is  called  a  phase  space.  Consider  a  two-dimensional  phase 
space  having  axes  labeled  by  the  components  of  the  dynamical  system  X  =  (Ar1,A'2)-  A 
point  in  the  phase  space  x  =  (xt  ,x2)  gives  a  complete  characterization  of  the  dynamical 
system  at  a  point  in  time.  As  time  proceeds  this  point  traces  out  a  curve  as  shown  in  Fig¬ 
ure  (2.1.1),  starting  from  the  initial  state  [X^OXX^O)]  and  proceeding  to  the  final  state 
[X  !(r),X2(r)]  at  time  t .  A  trajectory  or  orbit  in  phase  space  traces  out  the  evolution  of 
the  dynamical  system.  Time  is  a  parameter  which  indexes  each  point  along  such  a  solu¬ 
tion  curve.  The  field  of  trajectories  initiated  from  a  set  of  initial  conditions  is  often 
referred  to  as  the  flow  field.  If  for  example  the  flow  field  asymptotically  (r  — » °°)  con¬ 
verges  to  a  single  point  in  phase  space,  this  is  called  a  fixed  point  (or  focus)  [cf.  Fig¬ 
ure  (2.1.2)].  If  the  flow  field  converges  to  a  single  closed  curve  this  is  called  a  limit  cycle 
[cf.  Figure  (2.1.3)].  Such  limit  cycles  appear  as  periodic  behavior  in  the  variables  of 
interest. 

A  nonlinear  oscillator  which  is  “weakly”  nonlinear  is  capable  of  oscillating  at 
essentially  a  single  frequency  and  can  produce  a  signal  that  is  very  low  in  harmonic  con¬ 
tent.  Although  the  output  from  such  an  oscillator  system  is  sinusoidal  at  a  single  fre¬ 
quency,  there  are  fundamental  and  crucial  differences  between  such  an  oscillator  and  the 
classical  harmonic  oscillator,  the  latter  being  a  conservative  system  which  is  loss-free. 
The  basic  difference  is  that  the  nonlinear  oscillator  can  oscillate  at  one  and  only  one  fre¬ 
quency  and  one  and  only  one  amplitude,  the  amplitude  and  frequency  being  dependent 
on  one  another  for  a  given  configuration  of  parameters.  In  contrast,  the  amplitude  and 
frequency  are  independent  in  the  classical  linear  oscillator,  which  can  oscillate  at  any 
arbitrary  level  for  a  given  set  of  parameter  values.  These  differences  are  illustrated  in  the 
description  of  the  limit  cycle.  The  phase  plane  of  a  Hamiltonian  (loss-free)  oscillator  is 
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depicted  in  Figure  (2.1.4)  together  with  the  limit  cycle  for  an  oscillator  with  nonlinear 
dissipation  [cf.  Figure  (2.1.3)].  Although  there  are  superficial  resemblances  between 
these  diagrams,  there  are,  in  fact,  fundamental  differences  between  these  two  physical 
systems.  While  the  linear  conservative  oscillator  can  be  described  by  an  infinite  family 
of  closed  ellipses,  as  shown  in  Figure  (2.1.4),  the  nonlinear  oscillator  approaches  a  single 
limit  cycle  as  seen  in  Figure  (2.1.3).  This  limit  cycle  is  reached  asymptotically  whether 
the  initial  conditions  correspond  to  an  infinitesimal  perturbation  near  the  origin  or  to  a 
finite  perturbation  far  beyond  the  limit  cycle.  In  either  case  the  phase  point  spirals  to  the 
limit  cycle,  which  is  a  stable  final  state.  On  the  other  hand,  the  conservative  linear  oscil¬ 
lator  does  not  display  this  “structural  stability.”  Any  perturbation  causes  it  to  leave  one 
ellipse  and  move  to  another,  i.e.  the  orbits  are  neutrally  stable. 

In  linear  systems  the  term  equilibrium  is,  usually  applied  in  connection  with  conser¬ 
vative  forces,  with  the  point  of  equilibrium  corresponding  to  the  vanishing  of  all  forces 
with  the  system  being  at  rest.  The  stability  of  such  an  equilibrium  state  is  then  defined  by 
the  behavior  of  the  system  when  it  is  subject  to  a  small  perturbation  i.e.,  a  small  displace¬ 
ment  away  from  the  equilibrium  state  in  phase  space.  Roughly  speaking,  the  terms  sta¬ 
bility  and  instability  indicate  that  after  the  perturbation  is  applied  the  system  returns  to 
the  equilibrium  state  (stable)  or  that  it  continues  to  move  away  from  it  (unstable)  or  that 
it  does  not  move  (neutral  stability). 

(a)  Strange  attractors  (deterministic  chaos) 

The  appellation  “strange  attractor”  was  given  to  those  attractors  on  which,  unlike 
the  system  just  discussed,  the  system  dynamics  are  aperiodic.  This  means  that  a  deter¬ 
ministic  equation  of  motion  gives  rise  to  a  trajectory  whose  corresponding  time  series 
nowhere  repeats  itself  over  time;  it  is  chaotic.  The  term  “chaotic”  refers  to  the  dynam¬ 
ics  of  the  attractor,  whereas  “strangeness”  refers  to  the  topology  of  the  attractor.  From 
the  point  of  view  of  classical  statistical  mechanics  the  idea  of  randomness  has  tradition¬ 
ally  been  associated  with  the  weak  interaction  of  an  observable  with  the  rest  of  the 
universe.  The  traditional  view  requires  there  to  be  many  (an  infinite  number)  degrees  of 
freedom  that  are  not  directly  observed,  but  whose  presence  is  manifest  through  fluctua¬ 
tions  in  the  physical  observations  [see  eg.,  West  and  Lindenberg,  1987].  More  recently  it 
has  been  learned  that  in  a  nonlinear  system  with  even  a  few  degrees  of  freedom  chaotic 
motion  can  be  observed  [Lorenz,  1963]. 
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What  we  present  in  this  subsection  are  some  of  the  recent  results  obtained  in  non¬ 
linear  dynamics  that  lead  to  chaos.  First  we  briefly  review  the  classical  work  of  Lorenz 
(1963)  on  a  deterministic  continuous  dissipative  system  with  three  variables.  The  phase 
space  orbit  for  the  solution  to  the  Lorenz  system  is  on  an  attractor,  but  of  a  kind  on  which 
the  solution  is  aperiodic  and  therefore  strange.  We  discuss  this  family  of  aperiodic  solu¬ 
tions  and  discover  that  chaos  lurks  in  a  phase  space  of  dimension  three.  Rossler  (1978) 
points  out  that  if  oscillation  is  the  typical  behavior  of  low-dimensionai  dynamical  sys¬ 
tems,  then  chaos,  in  the  same  way,  characterizes  three-dimensional  continuous  systems. 

The  modem  view  of  randomness  discussed  in  the  Introduction  can  be  traced  back  to 
Poincare,  but  the  recent  avalanche  of  interest  dates  from  the  attempts  of  Lorenz  to  under¬ 
stand  the  short  term  variability  of  weather  patterns  and  thereby  enhance  their  predictabil¬ 
ity  [see  eg..  Halloway  and  West,  1984],  His  approach  was  to  represent  a  forced  dissipa¬ 
tive  geophysical  hydrodynamic  flow  by  a  set  of  deterministic  nonlinear  differential  equa¬ 
tions  with  a  finite  number  of  degrees  of  freedom.  By  forcing  we  mean  that  the  environ¬ 
ment  provides  a  source  of  energy  for  the  flow  field,  which  in  this  case  is  a  source  of  heat 
at  the  bottom  of  the  atmosphere.  The  dissipation  in  this  flow  extracts  energy  from  the 
temperature  gradient  but  the  forcing  term  puts  energy  back  in.  For  the  particular  physical 
problem  Lorenz  was  investigating,  the  number  of  degrees  of  freedom  he  was  eventually 
able  to  use  was  three,  let’s  call  them  X ,  Y ,  and  Z.  In  the  now  standard  form  these  equa¬ 
tions  are 


dX  v  v 

—  =  -  gX  +  csY 

d  X 

(2.1.2) 

—  =  -XY  +  r  X  -  Y 

(2.1.3) 

dX 

—  =XY  -bZ 

(2.1.4) 

d  X 

where  a,  r  and  b  are  parameters.  The  solutions  to  the  Lorenz  model  can  be  identified 
with  trajectories  in  phase  space.  What  is  of  interest  here  are  the  properties  of  the  non¬ 
periodic  bounded  solutions  in  this  three  dimensional  phase  space.  A  bounded  solution  is 
one  that  remains  within  a  restricted  domain  of  phase  space  as  time  goes  to  infinity. 

The  phase  space  for  the  set  of  equations  (2.1.2)  -  (2.1.4),  is  three-dimensional  and 
the  solution  to  them  traces  out  a  curve  Yt{x,y,z)  given  by  the  locus  of  values  of 
X(r)  =  [X(t),  y(f),Z(f)]  [cf.  Figure  (2.1.5)].  We  can  associte  a  small  volume 
V0(t)=Xo(t)Y0{t)ZQ(t )  with  a  perturbation  of  the  trajectory  and  investigate  how  this 
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volume  of  phase  space  changes  with  time.  If  the  original  flow  is  confined  to  a  region  R 

then  the  rate  of  change  of  the  small  volume  with  time  dV0 /dt  must  be  balanced  by  the 

•  • 
flux  of  volume  J(r)  =  V0(r)X(r)  across  the  boundaries  of  R.  The  quantity  X(r)  in  the 

flux  J  represents  the  time  rate  of  change  of  the  dynamical  variables  in  the  absence  of  the 

perturbations,  i.e.,  the  unperturbed  flow  field  that  can  sweep  the  perturbation  out  of  the 

region  R .  The  balancing  condition  is  expressed  by  an  equation  of  continuity  and  in  the 

physics  literature  is  written 

V0(t)  +  V-J(t)  =  0  (2.1.15) 

or  in  detail 

— L_  —  V0(r)  =  3X  X  +dY  Y  +3Z  Z  (2.1.16) 

V0(t)  dt  0  '  x  Y  z 

where  d/dt  (  =  9,  +  x  •  V^)  is  the  so-called  convective  or  total  derivative  of  the  volume. 
Using  the  equations  of  motion  (2.1.2)  -  (2.1.4)  in  (2.1.6)  we  obtain 

»--(.  +  »+».  (2.1.7) 

Equation  (2.1.7)  is  interpreted  to  mean  that  as  an  observer  moves  along  with  an  element 
of  phase  space  volume  V0(r)  associated  with  the  flow  field,  the  volume  will  contact  at  a 
rate  6+0  +  1,  i.e.,  the  solution  to  (2.1.7)  is  V0(r)  =  =  0)  exp  [  -  (6  +o  +  l)r].  Hence 

the  volume  goes  to  zero  as  t  -+°°  at  a  rate  which  is  independent  of  the  solutions 
X(t),Y(t)  and  Z(r).  As  pointed  out  by  Lorenz  ,  this  does  not  mean  that  each  small 
volume  shrinks  to  a  point  in  phase  space;  it  may  simply  become  flattened  into  a  surface, 
one  with  a  fractional  dimension,  ie.  a  non-integer  dimension  between  two  and  three. 
Consequently  the  total  volume  of  the  region  initially  enclosed  by  the  surface  R  shrinks  to 
zero  at  the  same  rate,  resulting  in  all  trajectories  become  asymptotically  confined  to  a 
specific  subspace  having  zero  volume  and  a  fractal  dimension  [Ott,  1985], 

To  understand  the  relation  of  this  system  to  the  kind  of  dynamical  situation  we  were 
discussing  in  the  preceding  section  we  must  study  the  behavior  of  the  system  on  the  lim¬ 
iting  manifold  to  which  all  trajectories  will  be  ultimately  confined.  This  cannot  be  done 
analytically  because  of  the  nonlinear  nature  of  the  equations  of  motion  for  the  Lorenz 
model.  Therefore,  these  equations  are  integrated  numerically  on  a  computer  and  the 
resulting  solution  is  depicted  as  a  curve  in  phase  space  for  particular  values  of  the 


495.chapter.2 


9- 1 2-'88 


-  17- 


parameters  cr,  b  and  r.  The  technical  details  associated  with  the  mathematical  under¬ 
standing  of  these  solutions  is  available  in  the  literature,  see  e.g.,  Ott  (1985)  or  Eckmann 
and  Ruelle  (1985)  and  of  course  the  original  discussion  of  Lorenz  (1963). 

In  Figure  (2.1.6)  we  display  the  behavior  of  Y(t)  for  3000  time  units.  After  reach¬ 
ing  an  early  peak  at  t  =  35,  Y  (t)  relaxes  so  a  relatively  stable  value  at  t  =  85  which  per¬ 
sists,  subject  to  systematically  amplified  oscillations,  until  near  r =1650.  Beyond  this 
time  Y  (t )  becomes  pulse-like  and  appears  to  change  signs  at  apparently  random  inter¬ 
vals.  This  irregularity  is  not  just  in  the  spacing  between  maxima  but  also  in  the  sign  of 
the  adjacent  maxima,  i.e.,  the  irregular  occurrence  of  a  number  of  peaks  of  one  sign 
before  a  peak  of  the  opposite  sign  occurs. 

In  Figure  (2.1.7a)  the  solution  manifold  in  the  three  dimensional  phase  space  is 
shown  and  (2.1.7b)  projects  the  solution  manifold  onto  the  (z  ,y  )-plane  and  the  (x  ,y  )- 
plane.  The  trajectory  indicated  is  not  complete,  but  is  that  segment  traversed  in  the  time 
interval  t  =  1400  to  1900.  The  points  C  and  C'  are  the  fixed  points  of  the  equations,  i.e., 
the  values  of  x,y ,  and  z  for  which  X  =  Y  =Z - Oin  (2.1.2)-(2.1.4),  which  for  r>l  yield 
X  -  Y  =  ±[b  [r -\)]m ,  Z  -  r — 1.  These  two  views  of  the  trajectory  indicate  that  the 
erratic  behavior  apparent  in  the  Y (t)  plot  [cf.  Figure  (2.1.6)]  arises  from  the  orbit  spiral¬ 
ing  around  one  of  the  fixed  points  C  orC'  for  some  arbitrary  period  and  then  jumping  to 
the  vicinity  of  the  other  fixed  point,  spiraling  around  that  for  a  while  and  then  jumping 
back  to  the  other  and  on  and  on.  If  the  number  of  times  the  orbit  circled  C  and  C’  were 
recorded  and  ordered,  the  resulting  sequence  would  be  random.  Virtually  all  trajectories 
finally  end  up  on  this  highly  unstable  manifold. 

The  strange  attractor  depicted  in  Figure  (2.1.5)  is  not  the  only  solution  to  the  Lorenz 
system  of  equations.  This  solution  was  obtained  for  the  parameter  values 
a  =  10, b  =  8/3, r  -  28.  If  the  values  a  =  10  and  b  =  8/3  are  held  fixed  and  r  is  increased 
from  zero,  a  wide  range  of  attractors  and  subsequent  dynamic  behaviors  are  obtained. 
The  possible  flow  patterns  make  the  transition  from  stable  equilibria  independent  of  ini¬ 
tial  conditions,  to  chaotic  attractors  that  are  sensitively  dependent  on  initial  conditions,  to 
“chaotic  transients’’  [Yorke  and  Yorke  1979]  in  which,  for  certain  initial  conditions,  an 
apparently  chaotic  trajectory  emerges  and  asymptotically  decays  into  a  stable  equilibria. 
The  decay  time  is  a  sensitive  function  of  the  initial  state. 
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Lorenz,  in  examining  the  solution  to  his  equations,  deduced  that  the  trajectory  is 
apparently  confined  to  a  surface.  Ott  (1985)  comments  that  the  apparent  “surface”  must  I 

have  a  small  thickness,  and  it  is  inside  this  thickness  that  the  complicated  structure  of  the 
strange  attractor  is  embedded.  This  is  where  the  folding  discussed  in  the  Introduction 
actually  occurs.  If  one  were  to  pass  a  transverse  line  through  this  surface,  the  intersec¬ 
tion  of  the  line  with  the  surface  is  a  set  of  dimension  D  with  0 <D  <1 .  The  structure  of  ® 

the  attractor  is  therefore  fractal,  and  the  stretching  and  folding  of  the  trajectory  discussed 
earlier  is  a  geometric  property  of  the  attractor. 

The  erratic  behavior  in  the  time  series  depicted  in  Figure  (2.1.6)  is  also  apparent  in  I 

the  associated  spectrum.  The  spectrum  is  the  mean  square  value  of  the  Fourier  transform 
of  a  time  series,  i.e.,  the  Fourier  transform  of  the  correlation  function.  Consider  the  solu¬ 
tion  X  (t );  it  will  have  a  Fourier  transform  over  a  time  interval  T  defined  by 


and  a  power  spectral  density  (PSD) 

I  Xj(tD)|  2 

S„«o)  =  Hm  ~  •  (2.1.9) 

In  Figure  (2.1.8)  we  display  the  power  spectral  densities  (PSD)  (to)  and  Szz  (to)  as  c  al- 
culated  by  Farmer,  Crutchfield,  Froehling,  Packard  and  Shaw  (1980)  using  the  trajectory 
shown.  It  is  apparent  from  the  power  spectra  density  using  the  X  (t)  time  series  that  there 
is  no  dominant  periodic  x -component  to  the  dynamics  of  the  attractor,  although  lower 
frequencies  are  favored  over  higher  ones.  The  power  spectral  density  for  the  Z(t)  time 
series  has  a  much  flatter  spectrum  overall,  but  there  are  a  few  isolated  frequencies  at 
which  energy  is  concentrated.  This  energy  concentration  would  appear  as  a  strong 
periodic  component  in  the  time  trace  of  Z(t ).  From  this  one  would  conclude  thatX(t)  is 
non-periodic,  but  that  Z(r)  possesses  both  periodic  and  non-periodic  components.  In  fact 
from  the  linearity  of  the  Fourier  transform  (2.1.8)  we  would  say  that  Z(r)  is  a  superposi¬ 
tion  of  these  two  parts: 

Z(t)  =  Zp(r)  +  Z„(t )  -  (2.1.10) 

The  implication  of  (2. 1.10)  is  that  the  auto-correlation  function 


r/2 


Xj (co)  =  j  X(r)e~iox  ~~ 
-  r/2  Z71 


(2.1.8) 
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C„ (t)  =  lim  <Z(f)Z(r+x)>  (2.1.11) 

may  be  written  as  the  sum  of  a  nonperiodic  components  <Znp(t)Znp(t  +x)>  that  decays 
to  zero  at  x->  °°  and  a  periodic  component  <  ZP  (t  )Zp  (t  +x)>  that  does  not  decay. 

To  summarize:  we  have  here  a  new  kind  of  attractor  that  is  referred  to  as  “strange’  ’ 
whose  dynamics  are  “chaotic”  and  with  a  power  spectra  density  resulting  from  the  time 
series  of  the  trajectory  that  has  broadband  components.  Dynamical  systems  that  are 
periodic  or  quasi-periodic  have  a  PSD  composed  of  delta  functions,  i.e.,  very  narrow 
spectral  peaks;  non-periodic  systems  have  broad  spectra  with  no  dramatic  emphasis  of 
any  particular  frequency.  It  is  this  broad  band  character  of  the  PSD  that  is  currently  used 
to  identify  non-periodic  behavior  in  experimental  data. 

So  what  does  this  all  mean ?  In  part  what  it  means  is  that  the  dynamics  of  a  complex 
system  might  be  random  even  if  its  description  can  be  “isolated”  to  a  few  (three  or 
more)  degrees  of  freedom  that  interact  in  a  deterministic  but  nonlinear  way.  If  the  sys¬ 
tem  is  dissipative,  i.e.,  information  is  extracted  from  the  system  on  the  average,  but  the 
system  is  open  to  the  environment,  i.e.,  information  is  supplied  to  the  system  by  means  of 
boundary  conditions,  then  a  “strange  attractor”  is  not  only  a  possible  manifold  for  the 
solutions  to  the  dynamic  equations;  it,  or  something  like  it,  may  even  be  probable. 

We  show  subsequently  that  the  aperiodic  or  chaotic  behavior  of  an  attractor  is  a 
consequence  of  a  sensitivity  to  initial  conditions:  trajectories  that  are  initially  nearby 
exponentially  separate  as  they  evolve  forward  in  time  on  a  chaotic  attractor.  Thus  as 
Lorenz  observed:  microscopic  perturbations  (unobservable  changes  in  the  initial  state  of 
a  system)  are  amplified  to  affect  macroscopic  behavior.  This  property  is  quite  different 
from  the  qualitative  features  of  nonchaotic  attractors.  In  the  latter,  orbits  that  start  out 
near  one  another  remain  close  together  forever.  Thus  small  errors  or  perturbations 
remain  bounded  and  the  behavior  of  individual  trajectories  remain  predictable. 

Of  course  these  considerations  are  not  of  much  practical  value  unless  they  can  be 
implemented  in  the  determination  of  the  properties  of  a  real  data  set.  This  will  be  done 
subsequently.  The  rationale  for  their  application  was  also  developed  by  Lorenz  in  his 
seminal  work,  but  the  full  extent  of  its  importance  has  only  recently  begun  to  emerge,  see 
e.g.  Lanford  (1976).  He  (Lorenz)  observed  that  the  trajectory  leaves  the  spiral  centered 
at  C  say,  [see  Figure  (2.1.5)],  only  after  exceeding  some  critical  distance  from  the  center. 
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Further,  the  degree  to  which  this  critical  distance  is  exceeded  determines  the  point  at 
which  the  next  spiral,  i.e.,  that  centered  at  C\  is  entered  as  well  as  the  number  of  circuits 
executed  prior  to  making  the  transition  back  to  the  C  center  again.  Thus  he  concludes 
that  “some  single  feature  of  a  given  circuit  should  predict  the  same  feature  of  the  follow¬ 
ing  circuit.”  As  an  example  he  selected  the  maximum  value  of  the  Z(t)  variable  along 
the  trajectory  which  occurs  whenever  the  circuit  is  nearly  completed. 

In  Figure  (2.1.9)  the  abscissa  is  labeled  by  the  value  of  the  nlh  maxima  Z„  of  Z(r) 
and  the  ordinate  is  labeled  by  the  value  of  the  following  maximum  Zn+1.  It  is  clear  that 
the  points  generated  lie  along  a  curve  if  the  spaces  between  points  are  filled  in.  This  is 
shown  for  example  by  Shaw  (1981)  using  the  increased  computing  capacity  that  has 
developed  in  the  intervening  years.  The  computer  generated  function  clearly  prescribe  a 
two-  to-one  relation  between  Z„  and  Zn+1.  From  this  relation  one  could  formulate  an 
empirical  prediction  scheme  using  the  geometry  of  the  attractor  as  a  data  set  without  a 
knowledge  of  the  underlying  dynamic  equations.  In  the  next  section,  after  we  learn  about 
mappings,  we  will  see  how  this  is  done. 

A  second  example  of  a  dynamic  system  whose  solutions  lie  on  a  chaotic  attractor 
was  given  by  Rossler  (1976)  for  a  chemical  process.  He  has  in  fact  provided  over  half  a 
dozen  examples  of  such  attractors  [cf.  Rossler  (1978)  ],  but  we  will  not  discuss  all  of 
them  here.  It  is  useful  to  consider  his  motivation  for  constructing  such  a  variety  of 
chaotic  attractors.  In  large  part  it  was  to  understand  the  detailed  effects  of  the  stretching 
and  folding  operations  in  nonlinear  dynamical  systems.  These  operations  mix  the  orbits 
in  phase  space  in  the  same  way  a  baker  mixes  bread  by  kneading  it,  i.e.,  rolling  it  out  and 
folding  it  over.  Visualize  a  drop  of  red  food  coloring  placed  on  top  of  a  ball  of  dough. 
This  red  spot  represents  the  initially  nearby  trajectories  of  a  dynamic  system.  Now  as  the 
dough  is  rolled  out  for  the  first  time  the  red  spot  is  stretched  into  an  ellipse,  which  even¬ 
tually  is  folded  over.  After  a  sufficiently  long  time  the  red  blob  is  stretched  and  folded 
many  times,  resulting  in  a  slab  of  dough  with  alternating  layers  of  red  and  white. 
Crutchfield  et  al.  (1987)  point  out  that  after  20  such  operations  the  initial  blob  has  been 
stretched  to  more  than  a  million  times  its  original  length,  and  its  thickness  has  shrunk  to 
the  molecular  level.  The  red  dye  is  then  thoroughly  mixed  with  the  dough,  just  as  chaos 
thoroughly  mixes  the  trajectories  in  phase  space  on  the  attractor. 
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The  dynamic  equations  for  Rossler’s  (1976)  three  degree  of  freedom  system  is 


X  =  -  (Y  +Z) 

(2.1.12) 

Y  =X  +aY 

(2.1.13) 

Z  =b  +XZ  -cZ 

(2.1.14) 

where  a,  b  and  c  are  constants.  For  one  set  of  parameter  values,  Farmer  et  al.  (1980) 
referred  to  the  attractor  as  “the  funnel,”  the  obvious  reason  for  this  name  is  seen  in  Fig¬ 
ure  (2.1.10).  Another  set  of  parameter  values  yields  the  “simple  Rossler  attractor,”  [cf. 
Figure  (2.1.1  Id)].  Both  of  these  chaotic  attractors  have  one  positive  Lyapunov  exponent. 
As  we  mentioned  earlier,  a  Lyapunov  exponent  is  a  measure  of  the  rate  at  which  trajec¬ 
tories  separate  one  from  the  other  (cf.  Section  (2.0)].  A  negative  exponent  implies  the 
orbits  approach  a  common  fixed  point.  A  zero  exponent  means  the  orbits  maintain  their 
relative  positions;  they  are  on  a  stable  attractor.  Finally,  a  positive  exponent  implies  the 
orbits  exponentially  separate;  they  are  on  a  chaotic  attractor. 

Equations  (2.1.19)  -  (2.1.21)  is  one  of  the  simplest  sets  of  differential  equation 
models  possessing  a  chaotic  attractor.  Figure  (2. 1.1 1)  depicts  a  projection  of  the  attractor 
onto  the  (x  ,y  )-plane  for  four  different  values  of  the  parameter  c.  Notice  that  as  c  is 
increased  the  trajectory  changes  from  a  simple  limit  cycle  with  a  single  maximum  [Fig¬ 
ure  (2.1.1  la)],  to  one  with  two  maxima  [Figure  (2.1.1  lb)]  and  so  on  until  finally  the  orbit 
becomes  aperiodic  [Figure  (2.1.1  Id)].  Here  again,  as  with  the  Lorenz  attractor,  we  can 
relate  the  nlh  maximum  of  say  X  (t)  to  the  (n  +  l)'rt  maximum.  This  can  be  done  by  not¬ 
ing  the  intersection  of  the  trajectory  in  Figure  (2.1.11)  to  a  line  placed  transverse  to  the 
attractor.  In  this  way  we  obtain  the  plot  of  the  maximum  shown  in  Figure  (2.1.12),  the 
curve  yielding  the  functional  equation  x„+1  =f(xn)  which  is  a  difference  equation.  This 
figure  suggests  how  we  can  replace  a  continuous  model  by  one  which  is  discrete.  We 
shall  return  to  this  procedure  in  the  following  section. 

2.2  Nonlinear  Mappings 

Now  that  we  have  seen  the  brand  of  chaos  that  a  continuos  strange  attractor  gives 
we  examine  a  one-dimensional  noninvertible  nonlinear  map.  This  mapping  is  the 
discrete  analog  of  the  logistic  equation  and  leads  to  a  subharmonic  bifurcation  of  the 
solution  eventually  resulting  in  a  kind  of  chaos  that  is  distinct  from  that  generated  by  the 
Lorenz  attractor.  The  results  are  quite  general  for  maps  with  a  single  maxima.  A  third 
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kind  of  chaos  is  related  to  that  found  on  the  Lorenz  attractor  in  that  a  two-dimensional 
invertible  map  is  shown  to  have  an  attractor  which  is  strange,  i.e.,  it  is  the  discrete  analog 
of  the  Lorenz  attractor. 

One  of  the  fascinating  aspect  of  these  maps  is  that  they  appear  to  be  the  natural  way 
in  which  to  describe  the  time  development  of  systems  in  which  successive  generations 
are  quite  distinct.  The  result  of  the  mathematical  analysis  is  that  for  certain  parameter 
regimes  there  are  a  large  number  of  classes  of  discrete  dynamical  models  (maps)  with 
chaotic  solutions.  The  chaos  associated  with  these  solutions  is  such  that  the  orbits  are 
periodic  or  erratic  in  time,  but  the  chaos  of  one  class  has  not  been  shown  to  be  the  same 
as  that  of  another  class.  However,  they  all  indicate  that  one  must  abandon  the  notion  that 
the  deterministic  nonlinear  evolution  of  a  process  implies  a  predictable  result.  One  may 
be  able  to  solve  the  discrete  equations  of  motion  only  to  find  a  chaotic  solution  that 
requires  a  distribution  function  for  making  predictions. 

Just  as  in  Section  2.1  we  wish  to  describe  the  dynamics  of  a  system  characterized  by 
an  N -component  vector  X  =  (XlrX2>  —  >XN)  and  again  in  order  to  determine  the  future 
evolution  of  the  system  from  its  present  state  we  must  specify  a  dynamic  rule  for  each  of 
the  components.  For  many  systems  the  variables  need  not  be  considered  continuous 
functions  of  time,  but  rather  to  be  functions  of  a  discrete  time  index  specifying  succes¬ 
sive  generations.  The  minimum  unit  of  time  change  for  the  dynamic  equations  would  in 
this  case  be  given  by  unity,  i.e.,  the  change  of  a  single  generation.  Thus  the  equations  of 
motion  instead  of  being  given  by  (2.1.1)  would  be  of  the  form 

X(n  +  1)  =  F  [  X(n )  ]  (2.2.1) 

where  the  changes  in  the  vector  X(n)  between  generation  n  and  n  f  1  is  determined  by 
the  function  F[X(n )].  If  at  generation  n  =0we  specify  the  components  of  X(0),  i.e.,  the 
set  of  circumstances  characterizing  the  system,  then  the  evolution  of  the  system  is  deter¬ 
mined  by  iteration  (mapping)  of  the  recursion  relation  (2.2.1.)  away  from  the  initial  state. 
Even  in  systems  that  are  perhaps  more  properly  described  by  continuous  time  equations 
of  motion  it  is  thought  by  many,  see  e.g.  Collete  and  Eckmann  (1980),  that  a  discrete 
time  representation  may  be  used  to  isolate  simplifying  features  a  certain  dynamical  sys¬ 
tems. 


495.chapter.2 


9-1 2-'88 


-23- 


(a)  One-dimensional  maps 

The  evolution  equation  in  a  discrete  representation  is  called  a  map  and  the  evolution 
is  given  by  iterating  the  map,  ie.,  by  repeated  application  of  the  mapping  operation  to  the 
newly  generated  points.  Thus  iterations  of  the  form  Xn  -»Xn+1  =/(X„),  where /  maps 
the  one-dimensional  interval  [0,1]  onto  itself,  is  interpreted  as  a  discrete  time  version  of  a 
continuous  dynamical  system.  The  choice  of  interval  [0,1]  is  arbitrary  since  the  change 
of  variables  Y  =  (X  -1  )/(b-a)  will  replace  a  mapping  of  the  interval  [a  ,b  ]  into  itself  by 
one  that  maps  [0,1]  into  itself.  For  example,  consider  the  continuous  trajectory  in  the 
two-dimensional  phase  space  depicted  in  Figure  (2.2.1).  The  intersection  points  of  the 
orb’t  with  the  X-axis  are  denoted  by  X1,X2,  •  •  •  .  The  point  Xn+1  can  certainly  be 
related  to  X„  by  means  of  the  function  /  determined  by  the  trajectory.  Thus,  instead  of 
solving  the  continuous  differential  equations  that  describe  the  trajectory,  in  this  approach 
one  produces  models  of  the  mapping  function  /  and  studies  the  properties  of 
X„+i=/(X„).  Here,  as  we  have  said,  n  plays  the  role  of  the  time  variables.  This  stra¬ 
tegy  has  been  applied  to  models  for  biological,  social,  economic,  chemical  and  physical 
systems.  May  (1976)  has  pointed  out  a  number  of  possible  applications  of  the  fundamen¬ 
tal  equation  for  a  single  variable 

X»+l  =/(*„)  •  (2.2.2) 

In  genetics,  for  example,  X„  could  describe  the  change  in  the  gene  frequency  between 
successive  generations;  in  epidemiology,  the  variable  X„  could  denote  the  fraction  of  the 
population  infected  at  time  n ;  in  psychology,  certain  learning  theories  can  be  cast  in  the 
form  where  X„  is  interpreted  as  the  number  of  bits  of  information  that  can  be  remem¬ 
bered  up  to  generation  n ;  is  sociology,  X„  might  be  interpreted  as  the  number  of  people 
having  heard  a  rumor  at  time  n  and  (2.2.2)  would  then  describe  the  propagation  of 
rumors  in  societies  of  various  structures  see,  e.g.,  Kemeny  and  Snell  (1972).  The  poten¬ 
tial  applications  of  such  modeling  equations  are  therefore  restricted  only  by  our  imagina¬ 
tions. 

Consider  the  simplest  mapping,  also  called  a  recursion  relation,  in  which  a  popula¬ 
tion  Xn  of  organisms  per  unit  area,  on  a  petri  dish  for  example,  in  the  nlh  generation  is 
strictly  proportional  to  the  population  in  the  preceding  generation  with  a  proportionality 
constant  ji; 


« 
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=M„_,  ,  n  =1,2,  •  (2.2.3) 

The  proportionality  constant  is  given  by  the  difference  between  the  birth  rate  and  death 
rate  and  is  therefore  the  net  birth  rate  of  the  population.  Equation  (2.2.3)  is  quite  easy  to 
solve.  Suppose  that  the  population  has  a  level  X0=N0  at  the  initial  generation,  then  the 
recursion  relation  yields  the  sequence  of  relation 

X{  =  \iN0  ,  X2  =  [LXl  =  ii2N0,-  -  (2.2.4) 


so  that  in  general 


Xn  =  pn  N0 


(2.2.5) 


This  rather  simple  solution  already  exhibits  a  number  of  interesting  properties.  Firstly,  if 
the  net  birth  rate  p  is  less  than  unity,  then  we  can  write  p"  -e  ~  where  (3>0,  so  that  the 
population  decreases  exponentially  between  successive  generation  (note  (3  =  -  In  p). 
This  is  a  reflection  of  the  fact  that  with  pel,  the  population  of  organisms  fails  to  repro¬ 
duce  itself  from  generation  to  generation  and  therefore  it  exponentially  approaches 
extinction: 

lim  Xn  =  0  if  p  <  1  .  (2.2.6) 

n  — v*> 

On  the  other  hand  if  p  >  1 ,  then  we  can  write  p”  =  e  n  ^  where  P  (  =  In  p)>  0,  so  the  popula¬ 
tion  increases  exponentially  from  generation  to  generation.  This  is  a  reflection  of  the  fact 
that  with  p>l  the  population  has  an  excess  at  each  generation  resulting  in  a  population 
explosion.  This  is  the  Malthus’  exponential  population  growth: 

lim  Xn  =<x>  if  p>i  .  (2.2.7) 

n  — >  <*> 

The  only  value  of  p  for  which  the  population  does  not  have  these  extreme  tendencies  is 
p=  1,  when,  since  the  population  reproduces  itself  exactly  in  each  generation,  we  obtain 
the  unstable  situation: 

lim  Xn=N0  (2.2.8) 

n  — >=« 
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Of  course  this  simple  model  is  no  more  valid  than  the  continuous  growth  law  of 
Malthus  (1798),  which  he  used  to  describe  the  exponential  growth  of  human  populations. 
A  more  scientifically  oriented  investigator,  Verhulst  (1844)  ,  put  forth  a  theory  that 
somewhat  mediated  the  pessimistic  view  of  Malthus.  Verhulst  noted  that  the  growth  of 
real  populations  is  not  unbounded.  He  argued  that  such  factors  as  the  availability  of 
food,  shelter,  sanitary  conditions,  etc.  all  restrict  (or  at  least  influence)  the  growth  of 
populations.  He  included  these  effects  by  making  the  growth  rate  p  a  function  of  the 
population  level.  His  arguments  allows  us  to  generalize  the  discrete  model  to  include  the 
effects  of  limited  resources.  In  particular,  the  birthrate  is  assumed  to  decrease  with 
increasing  population  in  a  linear  way: 

p-»p(X„)  =  p[l -Xn/Q]  (2.2.9) 

Where  ©  is  the  saturation  level  of  the  population.  Thus  the  linear  recursion  relation 
(2.2.3)  is  replaced  with  the  nonlinear  discrete  logistic  equation  , 

Xn+x  =  [iXn[\-Xn/@]  .  (2.2.10) 

It  is  clear  that  when  Xn<<®  the  population  grows  exponentially  since  the  nonlinear 
term  is  negligible.  However  at  some  point  the  ratio  Xn/Q  is  going  to  be  of  the  order 
unity  and  the  rate  of  population  growth  will  be  retarded.  When  Xn  =  ©  there  are  no  more 
births  in  the  population.  Biologically  the  regime  Xn>®  corresponds  to  a  negative 
birthrate,  but  this  does  not  make  biological  sense  and  so  we  restrict  the  region  of 
interpretation  of  this  model  to  [1  - Xn/® ]  >0.  Finally,  we  reduce  the  number  of  parame¬ 
ters  from  two,  (J.  and  0,  to  one  by  introducing  Yn  =Xn/@  the  fraction  of  the  saturation 
level  achieved  by  the  population.  In  terms  of  this  ratio  variable  the  recursion  relation 
(2.2.10)  becomes 

Yn+l=\iYn[\-Y„]  .  (2.2.11) 

Segal  (1984)  challenges  the  readers  of  his  book  (at  this  point  in  the  analysis  of  this  map¬ 
ping)  to  attempt  and  predict  the  type  of  behavior  manifest  by  the  solution  to  (2.2.1 1),  e.g. 
Are  there  periodic  components  to  the  solution?  Does  extinction  ever  occur?,  etc.  His 
intent  '''as  to  alert  the  reader  to  the  inherent  complexity  contained  in  the  deceptively  sim¬ 
ple  looking  equation  (2.2.1 1).  We  will  examine  some  of  these  general  properties  shortly, 
but  first  let  us  explore  our  example  a  bit  more  fully.  Our  intent  is  to  introduce  the  reader 
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to  a  number  of  fundamental  dynamical  concepts  that  will  be  useful  in  the  subsequent 
study  of  water  wave  data.  Note  that  saturation  inducing  models  such  as  (2.2.11)  have 
been  used  to  describe  the  transport  of  energy  in  physical  oceanography. 

We  noticed  that  extinction  was  the  solution  to  the  simple  system  (2.2.3)  when  p<l. 
Is  extinction  a  possible  solution  to  (2.2.11)?  If  it  is,  then  once  that  state  is  attained,  it 
must  remain  unchanged  throughout  the  remaining  generations.  Put  differently,  extinction 
must  be  a  steady-state  solution  of  the  recursion  relation.  A  steady-state  solution  is  one 
for  which  Yn  =Yn  +  l  for  all  n .  Let  us  assume  the  existence  of  a  steady-state  level  Yss  of 
the  population  such  that  (2.2.1 1)  becomes 

Yss=nYss(l-Yss)  (2.2.12) 

for  all  n,  since  in  the  steady-state  Yn+l  =  Yn  =  YSS.  Equation  (2.2.12)  defines  the  qua¬ 
dratic  equation 

yJ+(l/H- 1)^=0  (2.2.13) 

which  has  the  two  roots  Yss  =0,  and  Yss  =(1  -  l/|i).  The  Yss  =0  root  corresponds  to 
extinction,  but  we  now  have  a  second  steady  solution  to  the  mapping,  that  being 
Yss  =  l-  1/p.  One  of  the  questions  that  is  of  interest  in  the  more  general  treatment  of  this 
problem  is  to  determine  to  which  of  these  steady  states  the  population  evolves  as  time  go 
by,  ie.,  extinction  or  some  finite  constant  level. 

Before  we  examine  the  more  general  properties  of  (2.2.1 1)  and  equations  like  it,  let 
us  use  a  more  traditional  tool  of  analysis  and  examine  the  stability  of  the  two  steady 
states  found  above.  Traditionally  the  stability  of  a  system  in  the  vicinity  of  a  given 
value  is  determined  by  perturbation  theory.  We  use  that  technique  now  and  write 

(2-2.14) 

where  <  <  1  so  that  (2.2.14)  denotes  a  small  change  in  the  relative  population  from  its 
steady-state  value.  If  we  now  substitute  (2.2.14)  into  (2.2.1 1)  we  obtain 

Yss  +  Sn  +  l  =  \l(  Yss  +  )  [  1  -  Yss -U  •  (2.2. 15) 

Then  using  (2.2.12)  to  eliminate  certain  terms  and  neglecting  terms  quadratic  in  £,n  we 
obtain 

Ui  =  (2.2.16) 
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as  the  recursion  relation  for  the  perturbation.  In  the  neighborhood  of  extinction,  the 
*  Yss  =0  steady  state,  (2.2.16)  reduces  to  (2.2.3)  in  the  variable  rather  than  Xn.  There¬ 

fore  if  0  <p<l  then  the  fixed  point  Yss  =0  is  stable  and  if  p>l  the  fixed  point  is  unstable. 
By  stable  we  mean  that  % „  as  n  — ><»  if  0  <  pel  so  that  the  system  returns  to  the 

(  fixed  point,  i.e.,  decreases  exponentially  in  n .  By  unstable  we  mean  that  — »<»  as 

n  — *»  if  p>l  so  that  the  perturbation  grows  without  bound  and  never  returns  to  the  fixed 
point,  i.e.,  increases  exponentially  with  n.  Of  course  p  =  1  means  the  fixed  point  is 

neutrally  stable,  i.e.,  it  neither  return  to  nor  diverges  from  Yss  =  0. 
i 

In  the  neighborhood  of  the  steady  state  Yss  =  1  -  1/p.  the  recursion  relation  becomes 

Ut  =  (2-P&,  .  (2.2.17) 

The  preceding  analysis  can  again  be  repeated  with  the  result  that  if  1  >  2  -  p>  -  1  the 
fixed  point  Yss  =  1  -  1/p  is  stable  and  implies  that  the  birthrate  is  in  the  interval  l<p<3. 
The  stability  is  monotonic  for  l<p<2,  but  because  of  the  changes  in  sign  it  is  oscillatory 
for  2<p<3.  Similarly  the  fixed  point  is  unstable  for  0<p<l  (monotonic)  and  p>3  (oscilla- 
>  tory). 

Following  Olsen  and  Degn  (1985)  we  examine  the  nature  of  the  solutions  to  (2.2.1 1) 
as  a  function  of  the  parameter  p  a  bit  more  closely.  This  can  be  done  using  a  simple 

computer  code  to  evaluate  the  iterates  Yn.  For  0<p<4  insert  an  initial  value  0<T0<1 

I 

into  (2.2.1 1)  and  generate  a  Y x,  which  is  also  in  the  interval  10,1].  This  second  value  of 
the  iterate  is  then  inserted  back  into  (2.2.1 1)  and  a  third  value  Y  2  is  generated;  here  again 
0<T2<1.  This  process  of  generation  and  reinsertion  constitutes  the  dynamic  process, 
which  is  a  mapping  of  the  unit  interval  into  itself  in  a  two-  to-one  manner,  i.e.,  two 
values  of  the  iterate  at  step  n  can  be  used  to  generate  a  particular  value  of  the  iterate  at 
step  n+1.  In  Figure  (2.2.2a)  we  show  Yn  as  a  function  of  n  for  p  =  2.8  and  observe  that 
as  n  becomes  large  (n>10)  the  value  of  Yn  becomes  constant.  This  value  is  a  fixed  point 
of  the  mapping  equal  to  l-l/p  =  0.643,  and  is  approached  by  all  initial  conditions 
0<  K0<  1  i.e.,  it  is  an  attractor.  Quite  a  different  behavior  is  observed  for  the  same  initial 
point  when  p  =  3.2.  In  Figure  (2.2.2b)  we  see  that  after  an  initial  transient  the  process 
becomes  periodic,  that  is  to  say  that  the  iterate  alternates  between  two  values.  This 
periodic  orbit  is  called  a  2-cycle.  Thus,  the  fixed  point  becomes  unstable  at  the  parame¬ 
ter  value  p  =  3  and  bifurcates  into  a  2-cycle.  Here  the  2-cycle  becomes  the  attractor  for 
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the  mapping.  For  a  slightly  larger  value  of  }i,  jj. = 3.53,  the  mapping  settles  down  into  a 
pattern  in  which  the  value  of  the  iterate  alternates  between  two  large  values  and  two 
small  values  [cf.  Figure  (2.2.2c)].  Here  again  the  existing  orbit,  a  2-cycle,  has  become 
unstable  at  p.  =  3.444  and  bifurcated  into  a  4-cycle.  Thus,  we  see  that  as  (i  is  increased  a 
fixed  point  changes  into  a  2-cycle,  a  2-cycle  changes  into  a  4-cycle,  which  in  turn  will 
change  into  an  8-cycle  and  so  on.  This  process  of  period  doubling  is  called  subharmonic 
bifurcation  since  a  cycle  of  a  given  frequency  0)0  bifurcates  into  periodic  orbits  which 
are  subharmonics  of  the  original  orbit,  i.e.,  for  k  bifurcations  the  frequency  of  the  orbit  is 
coo/2*.  The  attractor  for  the  dynamic  process  can  therefore  be  characterized  by  the 
appropriate  values  of  ji. 

As  one  might  have  anticipated,  the  end  point  of  this  period  doubling  process  is  an 
orbit  with  an  infinite  period  (zero  frequency).  An  infinite  period  implies  that  the  system 
is  aperiodic,  that  is  to  say,  the  pattern  of  the  values  of  the  iterate  does  not  repeat  itself  in 
any  finite  number  of  iterations,  i.e.,  finite  time  interval,  [cf.  Figures  (2.2.2d)].  We  will 
see  presently  about  any  process  that  does  not  repeat  itself  as  time  goes  to  infinity  is  com¬ 
pletely  unique  and  hence  is  random.  It  was  this  similarity  of  the  mapping  to  discrete  ran¬ 
dom  sequences  that  motivated  the  coining  of  the  term  chaotic  to  describe  such  attractors. 
The  deterministic  mapping  (2.2.1 1)  can  therefore  generate  chaos  for  certain  values  of  the 
parameter  |i. 

Returning  now  to  the  more  general  context  it  may  appear  that  limiting  the  present 
analysis  to  one-dimensional  systems  is  unduly  restrictive;  however,  we  recall  that  the 
system  is  pictured  to  be  a  projection  of  a  more  complicated  dynamical  system  onto  a 
one-dimensional  subspace  [cf.  e.g..  Figure  (2.2.1)].  A  substantial  literature  based  on 
(2.2.11)  has  developed  in  the  past  decade,  much  of  which  is  focused  on  the  purely 
mathematical  properties  of  such  mappings.  The  physicists  and  mathematicians  have  been 
quite  actively  exploring  the  consequences  of  these  results  for  physical  and  chemical  sys¬ 
tems. 

For  the  moment  we  shall  make  the  assumption  that  the  maps  (dynamic  systems)  of 
interest  contain  a  single  maximum  and  that  f  (X)  is  monotonically  increasing  for  value  of 
X  below  this  maximum  and  monotonically  decreasing  for  values  of  X  above  this  max¬ 
imum.  Maps  such  as  these,  i.e.,  maps  with  a  single  maximum,  are  called  noninvertible, 
since,  given  Xn+1  there  are  two  possible  values  of  Xn  and  therefore  the  functional 
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relation  cannot  be  inverted.  If  the  index  n  is  interpreted  as  the  discrete  time  variable,  as 
we  did  above,  the  recursion  relation  generates  new  values  of  Xn  forward  in  time  but  not 
backward  in  time,  see  e.g.  Ott  (1985).  This  assumption  corresponds  to  the  reasonable 
requirement  that  the  dynamic  law  stimulates  X  to  grow  when  it  is  near  zero,  but  inhibits 
its  growth  when  it  reaches  some  saturation  value.  An  example  of  this  is  provided  by  the 
discrete  version  of  the  Verhulst  equation  for  population  growth  that  we  have  just  exam¬ 
ined.  Equation  (2.2.10)  is  often  called  the  discrete  logistic  equation  and  has  been  inten¬ 
sively  studied  in  the  physical  sciences,  usually  in  the  scaled  form  (2.2.11).  Thus  the 
mapping  function  is  f  (Yn)  =  \iYn(\  -Yn)  and  when  graphed  versus  Yn  yields  the  qua¬ 
dratic  curve  depicted  in  Figure  (2.2.3). 

The  mapping  operation  is  one  that  is  accomplished  by  applying  the  function  /  to  a 
given  initial  values  Y0  to  generate  the  next  point,  and  applied  sequentially  to  generate  the 
successive  images  of  this  point.  The  point  Yn  is  generated  by  applying  the  mapping  / ,  n 
times  to  the  initial  point  Y 0: 

Yn=fn(Y0)  (2.2.18) 

using  the  relation /"( T0)=/ {/n-1( T0)].  This  is  done  graphically  in  Figure  (2.2.3a)  for 
n  =3  using  the  rule:  starting  from  the  initial  point  T0  a  line  is  drawn  to  the  function 
yielding  the  value  Y  x  -f  ( T0)  along  the  ordinate,  then  from  symmetry  the  same  value  is 
obtained  along  the  abscissa  by  drawing  a  line  to  the  diagonal  (45°)  line.  An  application 
of/  to  T,  is  then  equivalent  to  dropping  a  line  from  the  diagonal  to  the  /  -curve  to  yield 
Y 2= f  ( T1)=/[/(To)]=/2(F0 ).  The  value  K3  is  obtained  in  exactly  the  same  way  from 
y3=/3(y0).  The  intersection  of  the  diagonal  with  the  function  /  defines  a  point  Y* 
having  the  property 

Y*  =  /  (T*)  (2.2.19) 

which  is  called  a  fixed  point  of  the  dynamic  equation,  i.e.,  Y*  is  the  Yss  from  (2.2.12). 
The  fixed  point  corresponds  to  the  steady-state  solution  of  the  discrete  equation  and  for 
(2.2.11)  Y*  =  1  -  1/p  (nontrivial)  and  Y*  =0  (trivial).  We  can  see  in  Figure  (2.2.3b)  that 
the  iterated  points  are  approaching  T*  and  as  n  — they  will  reach  this  fixed  point.  To 
determine  if  a  mapping  will  approach  a  fixed  point  asymptotically,  i.e.,  if  the  fixed  point 
is  stable,  we  examine  the  slope  of  the  function  at  the  fixed  point,  see  e.g..  May  (1976)  ; 
Li  and  Yorke  (1975)  ;  Collet  and  Eckmann  (1980).  The  function  acts  like  a  curved 
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mirror  either  focusing  the  ray  towards  the  fixed  point  under  multiple  reflections  or 
diverging  the  ray  away.  The  asymptotic  direction  (either  towards  or  away  from  the  fixed 
point)  is  determined  by  the  slope  of  the  function  at  Y* ,  which  is  depicted  in  Figure 
(2.2.4)  by  the  dashed  line  and  denoted  by  f'(Y*)  i.e.,  the  (tangent)  derivative  of  f  (Y)  at 
Y  =  Y* .  As  long  as  |  /'( Y*  )|  <1  the  iterations  of  the  map  are  attracted  to  Y  =Y* ,  just  as 
ihe  perturbation  2;  approaches  zero  in  (2.2.14)  near  the  stable  fixed  point.  Again  using  a 
logistic  map  as  an  example,  we  have  f'(Y*)=2-\i,  so  that  the  equilibrium  point  is 
stable  and  attracts  all  trajectories  originating  in  the  interval  0<F  <1  if  and  only  if  l<|i<3. 
This  is  of  course  the  same  result  we  obtained  using  linear  stability  theory  [cf.  Eq. 
(2.2.17)]  for  the  logistic  map. 

When  the  slope  of  f  is  such  that  the  fixed  point  becomes  unstable,  i.e.,  when 
\f\Y*)\  >1,  then  the  solution  “spirals”  out.  If  the  parameter  p  is  continuously 
increased  until  this  instability  is  reached  then  the  orbit  will  spiral  out  until  it  encounters  a 
situation  where  Y\  =/ (F*)  and  Y\=f  (Y \),  i.e.,  the  orbit  becomes  periodic.  Said  dif¬ 
ferently,  the  mapping  /  has  a  periodic  orbit  of  period  2  since  Y\  =f  (T* )  =/2(F 2 )  and 
F*  =/  (y^ )  =/2(F  * )  since  F*  and  Y\  are  fixed  points  of  the  mapping  / 2  and  not  of  the 
mapping  / .  In  Figure  (2.2.4a)  we  illustrate  the  mapping  /2  and  observe  it  to  have  two 
maxima  rather  than  the  single  one  of  / .  As  the  parameter  p  is  increased  further  the  dim¬ 
ple  between  the  two  maxima  increases  as  do  the  height  of  the  pealcs  along  with  the  slopes 
of  the  intersection  of/2  with  the  diagonal  [cf.  Figure  (2.2.4b)]. 

For  1<jj.<3  the  fixed  point  is  stable  and  Y*  is  a  degenerate  fixed  point  of  f2,  i.e., 
F*  =/2(F*).  At  p.  =  3.414  the  fixed  point  becomes  unstable  and  two  new  solutions  to 
the  quadratic  mapping  emerge.  These  are  the  two  intersections  of  the  quadratic  map  with 
the  diagonal  having  slopes  with  magnitude  less  than  unity,  F*  and  Y The  chain  rule  of 
differentiation  of  the  derivative  of  f2  at  Y\  and  Y 2  is  the  product  of  the  derivatives 
along  the  periodic  orbit 

f2\Y\ )=/'[/  (Y\ ))  f  \Y\ )  =  f  \Y\  )f  '(Y  2 )  =  f2'(Y  2 )  (2.2.20) 

so  that  the  slope  is  the  same  as  both  points  of  the  period  2  orbit,  see  e.g.,  Li  and  Yorke 
(1975) ,  and  in  fact  the  slope  is  the  same  at  all  k  of  the  values  of  a  period  k  orbit.  This  is 
in  fact  a  continuous  process  starting  from  the  stable  fixed  point  F*  when  |  f'\  <1;  as  p  is 
increased  this  point  becomes  unstable  at  |  / ']  =1  and  generates  two  new  stable  points 
with  \  f2'\  <1  for  a  period  2  orbit;  as  |i  is  increased  further  these  points  become  unstable 
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at  |  /  '|  =1  and  generates  two  new  stable  points  with  |  /2'|  <1  for  a  period  2  orbit;  as  p  is 
increased  further  these  points  become  unstable  at|  f2'\  =1  and  generates  four  new  stable 
points  with  |  f4'\  <1  for  a  period  4  orbit.  This  bifurcation  sequence  is  tied  to  the  value  of 
the  parameter  p.  As  this  parameter  is  increased  the  discrete  equation  undergoes  a 
sequence  of  bifurcations  from  the  fixed  point  to  stable  cycles  with  periods  2,  4,  8,  16,  32 
...  2*.  In  each  case  the  bifurcation  process  is  the  same  as  that  for  the  transition  from  the 
stable  fixed  point  to  the  stable  period  2  orbit.  A  graph  indicating  the  location  of  the 
stable  values  of  Y  for  a  given  p  is  given  in  Figure  (2.2.5).  Here  we  see  that  the  p  interval 
between  successive  bifurcations  is  diminishing  so  that  the  “window”  of  values  of  p 
wherein  any  one  cycle  is  stable  progressively  diminishes.  If  we  denote  by  p*.  the  value 
of  p  where  the  orbit  bifurcates  from  length  2i_1  to  2k,  then 

p*  -  p*.! 

lim  - =  universal  constant  (2.2.21) 

t  U*+i  “  U* 

a  result  first  obtained  numerically  by  Feigenbaum  (1979).  This  result  indicates  that  a 
constant  p„  is  being  approached  by  this  sequence.  This  critical  parameter  value  is  a 
point  of  accumulation  of  a  period  2k  cycles.  For  (2.2. 1 1)  the  critical  value  of  this  param¬ 
eter  is  Po,, =3.5700.  The  numerical  value  of  p^  is  dependent  on  the  particular  map  con¬ 
sidered,  although  the  existence  of  an  accumulation  point  does  not,  and  more  importantly 
the  universal  constant  in  (2.2.21)  has  a  value  4.69210  ...  and  is  also  independent  of  the 
specific  choice  of  the  map. 

In  Figure  (2.2.5)  we  use  the  logarithm  of  p  as  the  abscissa  in  order  to  clearly  distin¬ 
guish  the  bifurcation  points.  In  Figure  (2.2.6)  we  replot  this  sequence  linearly  in  p.  In 
the  latter  figure  we  distinguish  from  left  to  right,  a  stable  fixed  point,  orbit  of  period  1 ;  a 
stable  orbit  of  period  2,  then  4,  8  and  then  a  haze  of  orbits  starting  along  the  in  p^,  then 
another  orbit  of  period  6  then  5,  and  3.  Collet  and  Eckmann  (1980)  comment:  “The 
astonishing  fact  about  this  arrangement  of  stable  periodic  orbits  is  its  independence  of  the 
particular  one-parameter  family  of  maps.”  The  haze  of  points  beyond  p„  consists  of  an 
infinite  number  of  fixed  points  with  different  periodicities,  along  with  an  infinite  number 
of  different  periodic  orbits.  In  addition  there  are  an  uncountable  number  of  aperiodic  tra¬ 
jectories  (bounded)  each  of  which  is  associated  with  a  different  initial  point  Y 0.  Two 
such  adjacent  initial  points  generate  orbits  that  become  arbitrarily  distant  with  iteration 
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number;  no  mater  how  long  the  time  series  generated  by  f(Y)  is  iterated,  the  two  pat¬ 
terns  never  repeat.  As  mentioned,  Li  and  Yorke  (1975)  have  applied  the  term  chaotic  to 
this  hazy  region  where  an  infinite  number  of  different  trajectories  can  occur. 

Thus  we  have  arrived  at  the  remarkable  fact  that  a  simple  discrete  deterministic 
equation  can  generate  trajectories  that  are  aperiodic.  In  particular  in  order  to  form  a 
one-dimensional  map  to  exhibit  chaotic  behavior,  it  must  noninvertible. 

(h)  Two-dimensional  maps 

In  the  above  discussion  we  defined  a  mapping  in  terms  of  a  projection  of  a  higher 
order  dynamic  system  onto  a  one-dimensional  line.  This  same  definition  can  be  applied 
for  the  intersection  of  the  trajectories  of  a  higher  order  dynamic  process  with  a  two- 
dimensional  plane.  In  Figure  (2.2.7)  a  sketch  of  a  trajectory  in  three  dimensions  is 
shown,  the  intersection  of  the  orbit  with  a  plane  defines  a  set  of  points  that  can  be 
obtained  by  means  of  the  two-dimensional  map: 

Xn+i=fi<XH,Yn)  ,  Yn^=f2(Xn,Yn)  .  (2.2.22) 

Here  we  follow  Ott  (1°S5)  and  consider  only  invertible  maps  where  (2.2.22)  can  be 
solved  uniquely  for  Xn  and  Yn  as  functions  of  Xn+1  and  Kn+1;  Xn  -g\(Xn.  j,  Yn+1)  and 
Yn  =g2(-^n+i»  *n+i)-  If  n  is  the  time  index  then  invertibility  is  equivalent  to  time  rever¬ 
sibility,  so  that  these  maps  are  reversible  in  time  whereas  those  in  the  preceding  discus¬ 
sion  were  not.  The  maps  in  this  section  are  analogous  to  the  Hamiltonian  dynamic  equa¬ 
tions  discussed  in  physics  and  chemistry  and  not  the  dissipative  equations  leading  to  the 
strange  attractors  such  as  the  Lorenz  model. 

The  reason  for  examining  higher  order  maps,  such  as  the  two-dimensional  example 
given  by  (2.2.22)  is  that  under  certain  conditions  these  maps  have  many  of  the  properties 
of  the  so  called  strange  attractors  discussed  earlier  even  though  they  are  conservative.  In 
particular  we  will  establish  the  connection  between  these  invertible  maps  and  the  strange 
attractor  of  Lorenz  as  well  as  the  fractal  dimension  discussed  earlier.  > 

The  one -dimensional  noninvertible  maps  were  obtained  by  projecting  a  higher  order 
trajectory  onto  a  one-dimensional  line.  Let  us  now  reverse  the  process  and  expand  the 
space  of  the  noninvertible  map  from  one  to  two-dimensions  by  introducing  the  coordi- 
nate  Yn  in  the  following  way: 

*„  +  [  =f(Xn)  +  Yn  (2.2.23) 
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Of  course,  if  /  is  noninvertible  and  P=0  (2.2.23)  collapses  back  onto  the  one- 
dimensional  map  (2.2.11).  For  any  non-zero  p,  however,  the  map  (2.2.23)  is  invertible, 
i.e.,  Xn  =Fn+l/p  and  Yn  =Xn+i  -/ (T„+l/p).  Thus  we  have  transformed  a  noninvertible 
*  map  to  an  invertible  one  by  extending  the  space.  As  Ott  (1985)  points  out,  however,  if  P 

is  sufficiently  small  the  distinction  between  the  invertible  two-dimensional  map  and  the 
noninvertible  one-dimensional  map  may  not  be  measurable. 

Let  us  examine  the  behavior  of  a  small  phase  space  volume  as  the  two-dimensional 

> 

map  is  iterated  from  Xn  Yn  -  Vn  to  Xn+1  y„+1  =  Vn+1  in  analogy  to  what  was  done  with  the 
Lorenz  model.  The  relation  between  the  two  volumes  is 


i=JVn 

where  J  is  the  Jacobian  of  the  map: 

d^n+i  dXn+i 

ax,  ax, 

1  s  aiVw 

ax,  ax. 


(2.2.25) 


(2.2.26) 


Inserting  (2.2.23)  and  (2.2.24)  into  (2.2.26)  we  find  J  =  -  p  so  that  the  volume  at  con¬ 
secutive  times  (2.2.25)  is  given  by 


V*+\  =  -  PV, 


(2.2.27) 


which  for  an  initial  volume  V0  has  the  solution 

V„+1  =  (-l)n+1pn+1L0  ,  (2.2.28) 

so  that  if  |  p|  <  1  the  volume  will  contract  by  a  factor  [  PI  at  each  application  of  the  map¬ 
ping.  This  contraction  does  not  imply  that  the  solution  goes  over  to  a  point  in  phase 
space,  but  only  that  it  is  attracted  to  some  bounded  region  of  dimension  lower  than  that 
of  the  initial  phase  space.  If  the  dimension  of  the  attractor  is  non-integer,  then  the  attrac¬ 
tor  is  fractal;  see  e.g.  in  Mandelbrot  (1980)  where  the  observation  that  the  fractal  dimen¬ 
sion  of  a  set  may  or  may  not  be  consistent  with  the  term  strange.  Following  Eckmann 
(1981),  we  employ  the  property  that,  if  all  the  points  in  the  initial  volume  V0  converge  to 
a  single  attractor,  but  that  points  which  are  arbitrarily  close  initially  separate  exponen¬ 
tially  in  time,  then  that  attractor  is  called  strange.  This  property  of  nearby  trajectories  to 
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exponentially  separating  in  time  is  called  sensitive  dependence  on  initial  conditions  and 
gives  rise  to  the  aperiodic  behavior  of  strange  attractors.  There  exists  however  a  large 
variety  of  attractors  which  are  neither  periodic  orbits  nor  fixed  points  and  which  are  not 
strange  attractors.  All  of  these,  states  Eckmann  (1981),  seem  to  present  more  or  less  pro¬ 
nounced  chaotic  features.  Thus  there  are  attractors  that  are  erratic  but  not  strange.  We 
will  not  pursue  this  general  class  here. 

As  an  example  of  the  two-dimensional  invertible  mapping  we  first  transform  the 
logistic  equation  (2.2.1)  into  the  family  of  maps  Xn+1  =  1  -  cX2  with  the  parametric 
identification  c  ={yU2-  l)p/2  and  0<c  <2,  since  2<p.<4  and  Xn  maps  the  interval 
[-1,1]  onto  itself.  Then  using  (2.2.23)  and  (2.2.24)  we  obtain  the  mapping  first  studied 
by  Henon  (1976): 

Xn+l  =  l-cXn2  +  Yn  (2.2.29) 

Fn+1  =  PXB  (2.2.30) 

In  Figure  (2.2.8)  we  have  copied  the  loci  of  points  for  the  Henon  system  in  which  104 
successive  points  from  the  mapping  with  the  parameter  values  c  =  1.4  and  P  =  0.2  ini¬ 
tiated  from  a  variety  of  choice  of  (x0,y0).  Ott  (1985)  points  out  that,  as  the  map  is 
iterated,  points  come  closer  and  closer  to  the  attractor  eventually  becoming  indistinguish¬ 
able  from  it.  This,  however,  is  an  illusion  of  scale.  If  the  boxed  region  of  the  figure  is 
magnified  one  obtains  Figure  (2.2.9a)  from  which  a  great  deal  of  structure  of  the  attractor 
can  be  discerned.  If  the  boxed  region  in  this  latter  figure  is  magnified,  then  what  had 
appeared  as  three  unequally  space  lines  appear  in  Figure  (2.2.9b)  as  three  distinct  parallel 
intervals  containing  structure.  Notice  that  the  region  in  the  box  of  Figure  (2.2.9a) 
appears  the  same  as  that  in  Figure  (2.2.9b).  Magnifying  the  boxed  region  in  this  latter 
region  we  obtain  Figure  (2.2.9c),  which  aside  from  resolution  is  a  self-similar  representa¬ 
tion  of  the  structure  seen  on  the  two  preceding  scales.  Thus  we  observe  scale  invariant. 
Cantor-set-like  structure  transverse  to  the  linear  structure  of  the  attractor.  Ott  (1985) 
concludes  that  because  of  this  self-similar  structure  the  attractor  is  probably  strange.  In 
fact  it  has  been  verified  by  direct  calculation  that  initially  nearby  points  separate 
exponentially  in  time;  [see  eg.  Feit,  1978;  Curry,  1979],  thereby  coinciding  with  at  least 
one  definition  of  the  strange  attractor. 
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(c)  The  Lyapunov  exponent 

P  We  have  adopted  the  definition  that  chaotic  systems  are  those  that  have  a  sensitive 

dependence  on  initial  conditions.  This  sensitivity  requires  that  orbits  initially  near  to  one 
another  exponentially  separate  as  they  evolve  forward  in  time.  A  computable  quantita- 
^  tive  measure  of  the  rate  at  which  orbits  separate  is  the  Lyapunov  exponent.  For  a  one¬ 

dimensional  map  the  Lyapunov  exponent  is  defined  by  the  slope  of  the  map: 

1  N  ,  df(Yn) 

<7=  lim  jr  £  In  |  (2.2.31) 

n  — Pi  j  oi 

l 

where  T„+1=/  ( Yn ).  Shaw  (1981)  has  shown  that  c  is  also  the  average  information 
change  over  the  entire  interval  of  iteration.  He  argues  that  a  map  may  be  interpreted  as  a 
machine  that  takes  a  single  input  Y0  and  generates  a  string  of  numbers  during  the  itera- 
>  tion  process.  If  the  string  has  a  pattern  such  as  would  arise  for  an  attractor  that  is  a  fixed 

point  or  periodic  orbit,  then  after  a  very  short  time  the  machines  gives  no  new  informa¬ 
tion.  On  the  other  hand  if  the  orbit  is  chaotic  so  that  the  string  of  numbers  is  random, 
then  each  iterate  is  new  to  the  observer,  and  gives  a  new  piece  of  information.  Shaw 
convincingly  demonstrates  that  a  chaotic  process  is  a  generator  of  information.  He  argues 

that  a  negative  a  implies  a  periodic  orbit  and  the  magnitude  of  a  measures  the  degree  of 

stability  of  that  orbit  against  perturbations.  If  an  orbit  is  initiated  at  a  point  off  the 
t  periodic  orbit,  but  within  its  basin  of  attraction,  the  initial  data  will  be  lost  as  the  orbit 

damps  to  its  stable  values.  The  parameter  o  determines  the  rate  at  which  this  information 
is  lost  to  the  macroscopic  world.  If  a  is  positive,  then  it  determines  the  rate  of  diver¬ 
gence  of  nearby  trajectories  which  is  the  same  as  the  rate  of  information  production,  set 
1  Oseledec  (1968). 

As  an  example  let  us  take  p  =  4  in  (2.2. 1 1 ).  Then  if  we  define  a  new  variable 

Z„  =  —  sin-1  ( )  (2.2.32) 

K 

i 

the  logistic  map  transforms  to  the  “tent  map” 

f  2Zn  OSYn< 0.5 

Zn+l  = I  2(1  -Z„)  0.5<r„<l  '  (2.2.33) 

i  '■ 

From  this  we  obtain  for  the  slope  of  the  map  to  be  used  in  (2.2.31) 


> 
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for  all  Z ,  so  that 

a  =  ln2  «  0.693  >0  .  (2.2.35) 

Since  this  quantity  is  invariant  under  coordinate  transformations  this  proves  that  the 
logistic  map  with  p.  =  4  meets  our  definition  of  a  chaotic  dynamical  system,  i.e.,  0.693 
bits  of  information  are  generated  in  each  iteration.  In  fact  this  mapping  is  chaotic  for  all 
a  >  =  3.57  ■  •  •  since  o  >  0  for  all  these  values. 

Let  us  consider  an  N -dimensional  map,  i.e.,  X  =  (X  *,X2,  .  .  .  ,XN), 

X„+1=f(X„)  (2.2.36) 

for  which  we  have  a  trajectory  X'„  in  this  phase  space  with  initial  condition  X0  and  a 
nearly  trajectory  X„  with  initial  condition  X0  +  A  X0  and|  A  Xj  <  <|  X0|  .  Here  the  dou¬ 
ble  bars  denote  the  norm  of  the  vector.  The  difference  between  the  two  trajectories  A  X0 
defines  the  tangent  vector  u„  =  A  Xn  such  that  (2.2.36)  can  be  used  to  write 

u„+i  =  f(X'„)  -  f(X„)  (2.2.37) 

which  defines  the  linearized  mapping 

u„+i  =  A(Xn)  •  u„  .  (2.2.38) 

where  A  is  the  N  xN  matrix  defined  by 

A(Xn)  =  -|r  f(X„)  (2.2.39) 

so  that  the  map  (2.2.38)  is  linearized  along  the  trajectory  X„.  Following  Nicolis  (1986) 
the  solution  to  (2.2.38)  for  a  given  initial  condition  e0  at  the  nth  iteration  can  be  written 
as 

Un  =  Un,n-\  ^n-\,n-2  ^21^10^0  *  (2.2.40) 

where  U  is  the  fundamental  solution  matrix.  The  indexing  on  U  indicates  the  iteration 
for  which  it  is  the  solution  to  the  mapping.  Let  us  interpret  (2.2.40)  starting  with  the 
right-most  factor:  f/10e0  =  uI;  is  the  solution  (2.2.38)  for  the  initial  condition  X0.  The 
solution  u [  is  a  vector  of  length  d  (  and  director  et: 

U,0e0  =  d,e}  (2.2.41) 

and  ej  has  a  unit  norm.  Now  we  apply  f/2 1  to  e]  and  obtain  a  vector  of  length  d2  and 
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direction  e2.  Finally  we  can  rewrite  (2.2.40)  as  the  product  of  n  numbers 

un=d„d„_i  •••  dxen  ,  |e„|  =1  (2.2.42) 

instead  of  a  product  of  n  matrices.  The  maximal  Lyapunov  exponent  is  then  defined  as 

a  =  Iim  —  In  |  d  t  e„| 


n—*» o  n 


1  n 

=  lim  —  £  In  dn  . 
n—>  »=  n  •_  | 


(2.2.43) 


We  see  from  the  definition  of  dk  =J  Uk  jk_1  en_j|  that  In  dk  is  the  exponential  change  of 
the  length  of  e0  during  the  time  interval  when  the  system  (2.2.36)  moves  between  the 
iterates  X*_j  and  X*. 

Rather  than  finding  just  the  maximal  Lyapunov  exponent  we  can  define  a  Lyapunov 
exponent  for  each  of  the  N  variables  that  describe  the  dynamic  system.  To  do  this  we 
note  [cf.  Benettin,  Golgani  and  Strelcyn,  1976  ]  that  one  can  introduce  eigenvalues  \j(n ) 
of  the  matrix 


A„  =  [a(XJ  A(X„_!)  •••  A(X,)]!/n  , 


(2.2.44) 


where  A„  is  defined  by  (2.2.39)  and  is  the  Jacobian  matrix  of  f.  The  Lyapunov 
exponents  are  then  given  by 

a.  =  =  lim  In  |  X,-(n)|  (2.2.45) 

n  — >«*> 

These  eigenvalues  Xj  are  often  called  the  Lyapunov  numbers. 

Let  us  consider  the  example  given  by  Ott  (1985)  [cf.  Figure  (2.2.10)].  For  a  two- 
dimensional  map,  the  Lyapunov  numbers  are  given  by  Xj  and  X2  and  are  interpreted  as 
the  average  principle  stretching  factors  for  a  very  small  initial  circular  area  of  radius  e(0). 
More  formally  we  can  write 


\j  =  limj  magnitude  of  the  jth  eigenvaues  of 


r 

l/n  1 

A (Xn,Yn)  A(X„_„  Yn_0  • 

■  AIX^F,)]  > 

(2.2.46) 

where  A(X ,  Y)  is  the  Jacobian  matrix  of  the  map: 
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A(*,n 


a/i(x,n 

dx 

df2(X,Y) 

dX 


a/i(x,r) 

8K 

a/2(*.n 


ay 


(2.2.47) 


The  functions  / !  and  /  2  are  the  components  of  the  mapping  vector  f  in  (2.2.36);  and,  of 
course,  (Xt  Y  j),  •  •  •  ,  (Xn,  Yn)  is  a  sequence  generated  by  the  map.  Then  the  Lyapunov 
numbers  specify  the  average  stretching  rate  of  nearby  points.  If  the  map  is  to  be  chaotic, 
for  >A^  say,  then  must  be  greater  than  unity,  so  that  the  distance  between  almost 
nearby  points  increases  in  successive  iterations.  If  the  map  is  area  contracting  then 
Xi  X2  <1,  the  distance  between  almost  nearby  points  decreases  in  successive  iterations;  if 
it  is  area  preserving  then  ^^2=  1  and  the  distance  remains  unchanged. 


2.3  Measures  of  Strange  Attractors 

In  broad  outline  we  have  attempted  to  give  some  indications  of  how  simple  non¬ 
linear  dynamic  equations  can  give  rise  to  a  rich  variety  of  dynamic  behaviors.  In  particu¬ 
lar  we  have,  in  large  part,  focused  on  the  phenomenon  of  chaos  described  from  the  point 
of  view  of  mathematics  and  modeling,  but  little  or  no  effort  was  made  to  relate  these 
results  to  actual  data  sets.  Thus  the  techniques  may  not  appear  to  be  as  useful  as  they 
could  be  to  the  experimentalist  who  observes  large  variations  in  his/her  data  and  wonders 
if  the  observed  fluctuations  are  chaos  or  noise.  For  a  number  of  geophysical  phenomena 
there  may  be  no  reliable  dynamical  model  describing  the  behavior  of  the  system,  so  the 
investigator  must  use  the  data  directly  to  distinguish  between  the  two.  As  we  mentioned 
earlier,  a  traditional  method  for  determining  the  dynamic  content  of  a  time  series  is  to 
construct  the  power  spectrum  for  the  process  by  taking  the  Fourier  transform  of  the  auto¬ 
correlation  function,  or  equivalently  by  taking  the  Fourier  transform  of  the  time  series 
itself  and  forming  its  absolute  square  [cf.  (2.1.9)].  The  autocorrelation  function  provides 
a  way  to  use  the  data  at  one  time  to  determine  the  influence  of  the  process  on  itself  at  a 
latter  time.  It  is  a  measure  of  the  relation  of  the  value  of  a  random  process  at  one  instant 
of  time,  X  ( t )  say,  to  the  value  at  another  instant  x  seconds  later,  X{t  +  x).  If  we  have  a 
data  record  extending  continuously  over  the  time  interval  (-T  /2,  T 12),  then  the  autocorre¬ 
lation  function  is  defined  as 


Qx(T) 


lim 
T  — * 00 


T 


T/2 


J 

-77 2 


X  (t)X  ( t  +X)  dt 


(2.3.1) 
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Note  that  for  a  finite  sample  length,  i.e.  for  T  finite,  the  integral  defines  an  estimate  for 

the  autocorrelation  function  (x ,  T )  so  that  Ca (t)  =  lim  Ca( t,T).  In  Figure  (2.3.1) 

T  — 

a  sample  history  of  X(t)  is  given  along  with  its  displaced  time  trace  X  (t  +  x).  The  point 
by  point  product  of  these  two  series  is  given  in  (2.3.1)  and  then  the  average  over  the  time 
interval  ( -  772,  T H)  is  taken.  A  sine  wave,  or  any  other  harmonic  deterministic  data  set, 
would  have  an  autocorrelation  function  which  persists  over  all  time  displacements.  Thus 
the  autocorrelation  function  can  provide  a  measure  of  deterministic  data  embedded  in  a 
random  background. 

Similar  comments  apply  when  the  data  set  is  discrete  rather  than  continuous,  as  it 
would  be  for  the  mappings  in  Section  2.2.  In  the  discrete  case  we  denote  the  interval 
between  samples  as  A (  =  T IN)  for  N  equally  spaced  intervals  and  r  as  the  lag  or  delay 
number  so  that  the  estimated  autocorrelation  function  is 

C«(rA,N)  =  — i-  N£  X:  Xj+r  ,  r  =0,  1,  •  ■  •  ,  m  (2.3.2) 

and  m  is  the  maximum  lag  number.  Note  that  Cw(r  A  ,N)  is  analogous  to  the  estimate  of 
the  continuum  autocorrelation  function  and  becomes  the  true  autocorrelation  function  in 
the  limit  N—>° These  considerations  have  been  discussed  at  great  length  by  Wiener 
(1949)  in  his  classic  book  on  time  series  analysis,  and  is  still  recommended  today  as  a 
►  text  from  which  to  capture  a  master’s  style  of  investigation. 

The  frequency  content  is  extracted  from  the  autocorrelation  function  by  applying  a 
filter  in  the  form  of  a  Fourier  transform.  This  yield  the  power  spectral  density 

oo 

S„(a»  =  j  J  e  ~  i's>tCxx(t)dt  (2.3.3) 


of  the  time  series  X  (t).  Equation  (2.3.3)  relates  the  autocorrelation  function  to  the  power 
spectral  density  and  is  known  as  the  Weiner-Khinchine  relation  which  is  in  agreement 
with  (2.1.8).  One  example  of  its  use  is  provided  in  Figure  (2.3.2a)  where  the  exponential 
form  of  the  autocorrelation  function  Cxx(t)  =  e  ~I/Zc  used  in  Figure  (2.3.2b)  yields  a  fre¬ 
quency  spectrum  of  the  Cauchy  form 


*«(<*» =4- 


K  1  +  C02X2 


At  high  frequencies  the  spectrum  (2.3.4)  is  seen  to  fall-off  as  (o' 


(2.3.4) 
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As  we  mentioned,  a  periodic  signal  in  the  data  will  show  sharp  peaks  in  the  spec¬ 
trum  corresponding  to  the  fundamental  frequency  and  its  higher  harmonics  [cf.  Section 
2.1(b)].  On  the  other  hand  the  spectrum  corresponding  to  aperiodic  variations  in  the  time 
series  will  be  broadband  in  frequency  with  no  discernible  structure.  In  themselves  spec¬ 
tral  techniques  have  no  way  of  discriminating  between  chaos  and  noise  and  are  therefore 
of  little  value  in  determining  the  source  of  the  fluctuations  in  a  data  set.  They  were  in 
fact  very  useful,  as  shown  in  Section  2.1(b),  in  establishing  the  similarities  between  sto¬ 
chastic  processes  and  chaos  defined  as  the  sensitive  dependence  on  initial  conditions  in  a 
dynamic  process. 

One  way  in  which  some  investigators  have  proceeded  in  discriminating  chaos  from 
noise  is  to  visually  examine  time  series  for  period  doublings.  This  is  a  somewhat  risky 
business,  however,  and  may  lead  to  misinterpretations  of  data  sets.  Also,  period  doubling 
is  only  one  of  the  possible  routes  to  chaos  in  dynamic  systems. 

It  has  been  suggested  that  fractal  processes  associated  with  scaled,  broadband  spec¬ 
tra  are  “information-rich.”  Periodic  states,  in  contrast,  reflect  narrow-band  spectra  and 
are  defined  by  monotonous,  repetitive  sequences,  depleted  of  information  content.  In 
Figure  (2.3.3)  we  depict  the  spectrum  of  the  time  series  X  (r )  obtained  from  the  funnel 
attractor  solution  of  the  equation  set  (2.1.2)-(2.1.14).  The  attractor  itself  is  shown  in  Fig¬ 
ure  (2.1.10).  We  see  that  the  spectrum  is  broad  band  as  was  that  of  the  Lorenz  attractor 
[cf.  Figure  (2.1.8)],  with  a  number  of  relatively  sharp  spikes.  These  spikes  are  manifesta¬ 
tions  of  a  strong  periodic  components  in  the  dynamics  of  the  funnel  attractor.  Thus  the 
dynamics  could  easily  be  interpreted  in  terms  of  a  number  of  harmonic  components  in  a 
noisy  background,  but  this  would  be  an  error.  One  way  to  distinguish  between  these  two 
interpretations  is  by  means  of  the  information  dimension  of  the  time  series.  The  dimen¬ 
sion  decreases  as  a  system  undergoes  a  transition  from  chaotic  to  periodic  dynamics. 

Thus  we  conclude  that  more  systematic  methods  for  distinguishing  between  chaos 
and  noise  are  desirable  and  necessary.  We  turn  to  those  methods  now. 

(a)  Correlational  dimension 

In  the  preceding  discussion  we  presented  the  standard  example  of  a  correlation 
function  having  an  exponential  form.  Such  a  correlation  function  could  describe  a  ran¬ 
dom  time  series  having  a  memory  or  correlation  time  tc .  It  could  not  describe  a  dynami¬ 
cal  system  having  an  asymptotic  stationary  or  periodic  state.  Similarly  it  could  not 
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describe  a  nonlinear  dissipative  dynamical  system  that  has  a  chaotic  attractor. 
Grassberger  and  Procaccia  (1983a,  b,  c)  developed  a  correlational  technique  by  which 
one  can  exclude  various  choices  for  the  kind  of  attractor  on  which  the  dynamics  for  a 
given  data  set  exists.  They  wanted  to  be  able  to  say  that  the  attractor  for  the  data  set  is 
not  multiply  periodic,  or  that  the  irregularities  are  not  due  to  external  noise,  etc.  As  we 
have  just  seen  Fourier  analysis  would  tell  us  if  the  attractor  were  multiply  periodic,  but 
not  the  source  of  the  fluctuations.  They  proposed  a  measure  obtained  by  considering 
correlations  between  points  of  a  time  series  taken  from  a  trajectory  on  the  attractor  after 
the  initial  transients  have  died  away. 

Consider  the  set  (X,  ,j=  1,2,  ■  ■  • ,  N  }  of  points  on  the  attractor  taken  from  a  time 
scries  X  (/),  i.e.,  we  take  X;  =  X(jA)  where  A  is  a  fixed  time  interval  between  successive 
measurements.  We  see  that  this  set  of  points  could  also  be  determined  from  a  mapping 
where  j  denotes  the  iterate  of  the  map.  If  the  attractor  is  chaotic  then  since  nearby  trajec¬ 
tories  exponentially  separate  in  time,  we  expect  that  most  pairs  of  points  Xy ,  Xk  j*k  will 
be  dynamically  uncorrelated.  Even  though  these  points  may  appear  to  be  essentially  ran¬ 
dom,  they  do  all  lie  on  the  same  attractor  and  therefore  are  correlated  in  phase  space.  As 
we  discuss  in  the  next  section,  the  single  time  series  can  be  used  to  construct  an  m- 
dimensional  representation  of  the  data  by  shifting  the  data  as  follows: 
Xy(1)  =  X  (j A);Xy(2/=X  (y'A-x),  .  .  . ,  Xy(m)=X  [j A-  (m  -  l)x],  where  i  has  a  value 
that  is  after  determined  by  trial  and  error.  The  single  time  series  Xj  is  thus  replaced  by 
the  vector  time  series  X;  =  {Xj^,Xj®, . . . ,  Xfm^}  where  m  is  called  the  embedding 
dimension. 

Grassberger  and  Procaccia  (1983)  introduced  the  correlation  integral  C(r)  defined 
by 

c"(r)  *  Jim  ttt  X  e(r-|xj-x>l ) 

iV->~  N  ij-  1 
r 

=  jdmr'c(r')  (2.3.5) 

0 

where  0(x )  is  the  Heaviside  function,  =  0  if  x  SO  and  =  1  if  x  >  0,  and  c  ( r' )  is  the  tradi¬ 
tional  correlation  function  in  Euclidian  volume  of  m  -dimensions 

c(r)=  lim  — L  £  5m  (  X;  -  X.  -  r )  .  (2.3.6) 

N-*~N2 
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and  |  X,  -  \j\  is  the  Euclidean  norm  of  the  m  -dimensional  vector.  The  virtue  of  the 
integral  function  is  that  for  a  chaotic  or  strange  attractor  the  correlational  integral  has  the 
power- law  form 

lim  Cm(r)~rv  (2.3.7) 

m  — 

and  moreover,  the  “correlation  exponent”  v  is  closely  related  to  the  fractal  dimension  D 
and  the  information  dimension  a  of  the  attractor.  They  argue  that  the  correlation 
exponent  is  a  useful  measure  of  the  local  properties  of  the  attractor  whereas  the  fractal 
dimension  is  a  purely  geometric  measure  and  is  rather  insensitive  to  the  local  dynamic 
behavior  of  the  trajectories  on  the  attractor  The  information  dimension  is  somewhat  sen¬ 
sitive  to  the  local  behavior  of  the  trajectories  and  is  a  lower  bound  on  the  Hausdorff 
dimension.  In  fact  they  observe  that  in  general  one  has 

v  <o<D  .  (2.3.8) 

Thus  if  the  correlation  integral  obtained  from  an  experimental  data  set  has  the  power-law 
form  (2.3.7)  with  v  <m ,  one  knows  that  the  data  set  arises  from  deterministic  chaos 
rather  than  random  noise,  because  noise  will  result  in  Cm(r)~rm  for  a  constant  correla¬ 
tion  function  over  the  distance  r .  Note  that  for  periodic  sequences  v  =  1 ;  for  random 
sequences  it  should  equal  the  embedding  dimension,  while  for  chaotic  sequences  it  is 
finite  and  non-integer. 

Grassberger  and  Procaccia  (1983)  point  out  that  one  of  the  main  advantages  of  the 
correlation  dimension  v  is  the  ease  with  which  it  can  be  measured.  In  particular  it  can  be 
measured  more  easily  than  either  a  or  D  for  cases  when  the  fractal  dimension  is  large 
(>  3).  Just  as  they  anticipated,  the  measure  v  has  proven  to  be  most  useful  in  experimen¬ 
tal  situations,  where  typically  high  dimensional  systems  exist. 

To  test  their  ideas  they  studied  the  behavior  of  a  number  of  simple  models  for  which 
the  fractal  dimension  is  known.  In  Figure  (2.3.4)  we  display  three  of  the  many  calcula¬ 
tions  they  did.  In  each  case  the  logarithm  of  the  correlation  integral  is  plotted  as  a  func¬ 
tion  of  the  logarithm  of  a  dimensionless  length  which  according  to  the  power-law  rela¬ 
tion  (2.3.7)  should  yield  a  straight  line  of  positive  slope.  The  slope  of  the  line  is  the 
correlational  dimension  v.  We  see  from  these  examples  that  the  technique  successfully 
predicts  the  correlational  behavior  for  both  mappings  and  differential  equations  having 
chaotic  attractors. 
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(b)  Attractor  reconstruction  from  data 

More  often  than  not  the  geophysical  observer  does  not  have  the  luxury  ot  a 
mathematical  model  to  guide  the  measurement  process.  What  is  usually  available  are  a 
few  partial  theories,  securely  based  on  assumptions  often  made  more  for  convenience 
than  for  reality,  and  a  great  deal  of  phenomenology.  Therefore  in  a  system  known  to 
depend  on  a  number  of  independent  variables  it  is  not  clear  how  many  kinds  of  measure¬ 
ments  one  should  make.  In  fact  it  is  often  unrealistically  difficult  to  take  more  than  the 
measurement  of  a  single  degTee  of  freedom.  What  then  can  one  say  about  a  complex  sys¬ 
tem  given  this  single  time  series?  It  turns  out  that  quite  a  lot  may  be  learned  using 
methods  developed  in  nonlinear  dynamics.  In  particular  a  method  has  been  devised  that 
enables  one  to  reconstruct  a  multidimensional  attractor  from  the  time  series  of  a  single 
observable.  The  application  of  this  technique  to  a  number  of  geophysical  data  sets  will 
be  reviewed  in  the  next  chapter,  but  for  the  moment  we  concentrate  on  the  exposition  of 
the  underlying  theory. 

Packard,  Crutchfield,  Farmer  and  Shaw  (1980)  who  constituted  the  nucleus  of  the 
Dynamic  Systems  Collective  at  the  University  of  California,  Santa  Cruz  in  the  late  70’ s 
and  early  80’s,  were  the  first  investigators  to  demonstrate  how  one  reconstructs  a  chaotic 
attractor  from  an  actual  data  set.  They  used  the  time  series  generated  by  one  coordinate 
of  the  three-dimensional  chaotic  dynamical  system  studied  by  Rossler  (1978)  i.e., 
(2.1. 12)-(2. 1. 14)  with  the  parameter  values  a  =  0.2,  b  =  0.4  and  c  =  5.7.  The  reconstruc¬ 
tion  method  is  based  on  the  hueristic  idea  that  for  such  a  three-dimensional  system,  any 
three  “independent”  time  varying  quantities  are  sufficient  to  specify  the  state  of  the  sys¬ 
tem.  The  three  dynamic  coordinates  X(t),Y(t)  and  Z  (r)  are  only  one  of  the  many  possi¬ 
ble  choices.  They  conjectured  that;  “any  such  sets  of  three  independent  quantities  which 
uniquely  and  smoothly  label  the  states  of  the  attractor  are  diffeomorphically 
equivalent.”  In  English  this  means  that  an  actual  dynamic  system  does  not  know  of  the 
particular  representation  chosen  by  us,  and  that  any  other  representation  containing  the 
same  dynamic  information  is  just  as  good.  Thus,  an  experimentalist  sampling  the  values 
of  a  single  coordinate  need  not  find  the  “one”  representation  favored  by  nature,  since 
this  “one”  does  not  in  all  probability  exist. 
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Packard  et  al.  (1980)  playing  the  role  of  experimentalists  sampled  the  X(r)  coordi¬ 
nate  of  the  Possler  attractor.  They  then  noted  a  number  of  possible  '’Itematives  to  the 
phase  space  coordinates  (x,y,z)  that  could  give  a  faithful  representation  of  the  dynamics 
using  the  time  series  they  had  obtained.  One  possible  set  was  the  time  series  itself  plus 
two  replicas  of  it  displaced  in  time  by  x  and  2x,  i.e.  X  (r ),  X  (t -x)  and  X  (r  -2x).  Note  that 
implicit  in  this  choice  is  the  idea  that  X  (t )  is  so  strongly  coupled  to  the  other  degrees  of 
freedom  that  it  contains  dynamic  information  about  these  coordinates  as  well  as  itself.  A 
second  representation  set  could  be  obtained  by  making  the  time  interval  x  an 
infinitesimal,  so  that  by  taking  differences  between  the  variables  we  would  obtain 
X(t),X(t)  andX(r). 

Figure  (2.3.5a)  shows  a  projection  of  the  Rossler  chaotic  attractor  on  the  (x ,  y ) 
plane.  Figure  (2.3.5b)  depicts  the  reconstruction  of  that  attractor  from  the  sampled X(t) 
time  series  in  the  (x,x)  plane.  It  is  clear  that  the  two  attractors  are  not  identical,  but  it  is 
just  as  clear  that  the  reconstructed  one  retains  the  topological  characteristics  and 
geometrical  form  of  the  experimental  attractor.  One  quantitative  measure  of  the 
equivalence  of  the  experimental  and  reconstructed  attractors  is  the  Lyapunov  exponent 
associated  with  each  one.  This  exponent  can  be  determined  by  constructing  a  return  map 
for  each  of  the  attractors  and  then  applying  the  relation  (2.2.31). 

A  return  map  is  obtained  by  constructing  a  Poincard  surface  of  section.  In  this 
example  of  an  attractor  projected  onto  a  two-dimensional  plane,  the  Po incard  surface  of 
section  is  the  intersection  of  the  attractor  with  a  line  transverse  to  the  attractor.  We  indi¬ 
cate  this  by  the  dashed  line  in  Figure  (2.3.5b)  and  the  measured  data  are  the  sequence  of 
values  {X„  }  denoting  the  crossing  of  the  line  by  the  attractor  in  the  positive  direction. 
These  data  are  used  to  construct  a  next  amplitude  plot  in  which  each  amplitude  X„+1  is 
plotted  as  a  function  of  the  preceding  amplitude  Xn .  It  is  possible  for  such  a  plot  to  yield 
anything  from  a  random  spray  of  points  to  a  well  defined  curve.  If  in  fact  we  find  a  curve 
with  a  definite  structure  then  it  may  be  possible  to  construct  a  return  map  for  the  attrac¬ 
tor.  For  example,  the  oscillating  chemical  reaction  of  Belousov  and  Zhabotinskii  was 
shown  by  Simoyi,  Wolf,  and  Swinney  (1982)  to  be  describable  by  such  a  one- 
dimensior.al  map.  In  Figure  (2.3.6)  we  indicate  the  return  map  constructed  from  the 
experimental  data  of  Simoyi  et  al.  (1980),  also  Figure  (2.1.9)  for  the  Lorenz  attractor. 
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Simoyi  et  al.  (1982)  point  out  that  there  are  25  distinct  chemicals  in  the  Belousov- 
Znabotinskii  reaction,  many  more  than  can  be  reliably  monitored.  Theiefoie  there  is  no 
practical  way  to  construct  the  twenty-five  dimensional  phase  space 
m- 1,  ■  •  •  ,25  from  the  experimental  data.  Instead  they  use  the  embedding  theorems  of 
Whitney  (1936)  and  Takens  (1981)  to  justify  the  monitoring  of  a  single  chemical  species, 
in  this  case  the  concentration  of  the  bromide  ion,  for  use  in  constructing  an  m- 
dimensional  phase  portrait  of  the  attractor  (X(r),  X(t  -x),  •••  X[r-(m-l)x]}  for 
sufficiently  large  m  and  for  almost  any  time  delay  x.  They  find  that  for  their  experimen¬ 
tal  data  m=  3  is  adequate  and  the  resulting  one -dimensional  map  [cf.  Figure  (2.3.6)]  pro¬ 
vided  the  first  example  of  a  physical  system  with  many  degrees  of  freedom  that  can  be  so 
modeled  in  detail. 

Let  us  now  recap  the  technique.  We  assume  that  the  system  of  interest,  a  geophysi¬ 
cal  time  series  say,  can  be  described  by  m  variables,  where  m  is  large  but  unknown,  so 
that  at  any  instant  of  time  there  is  a  point  X(t)  =  (X^(r),  ■  •  •  ,X^m\r))  in  an 

m  -dimensional  phase  space  that  completely  characterizes  the  system.  This  point  moves 
around  as  the  system  evolves,  in  some  cases  approaching  a  fixed  point  or  limit  cycle 
asymptotically  in  time.  In  other  cases  the  motion  appears  to  be  purely  random  and  one 
must  distinguish  between  a  system  confined  to  a  chaotic  attractor  and  one  driven  by 
noise.  In  experiments,  one  often  only  records  the  output  of  a  single  detector,  which 
selects  one  of  the  m  components  of  the  system  for  monitoring.  In  general  the  experimen¬ 
talist  does  not  know  the  size  of  the  phase  space  since  the  important  dynamic  variables  are 
usually  not  known  and  therefore  he/she  must  extract  as  much  information  as  possible 
from  the  single  time  series  available,  X^(r)  say.  For  sufficiently  long  delay  times  x  one 
uses  the  embedding  theorem  to  construct  the  sequence  of  displaced  time  series 
{X(1)(r),X(1)(r  +x),  •  •  •  ,  X(1)  [t  -(m  -  l)x]  }  .  This  set  of  variables  has  been  shown  to 
have  the  same  amount  of  information  as  the  m  -dimensional  phase  point  Thus,  as  time 
goes  to  infinity,  we  can  build  from  the  experimental  data  a  one-dimensional  phase  space 
X^(f ),  a  two-dimensional  phase  space  with  axes  {X(t)(r), X(1\r -x)},  and  so  on.  One 
then  determines  if  the  system  dynamics  are  confined  to  a  low-dimensional  attractor. 
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3.  Review  of  Some  Geophysical  Applications  of  the  Reconstruction  technique  and 
other  Nonlinear  Concepts 

Ip  the  previous  Sections  we  have  made  every  effort  to  develop  some  rather  difficult 
mathematical  concepts  and  techniques  in  a  context  that  would  make  their  importance 
self-evident  in  a  geophysical  setting.  In  the  present  section  we  review  how  a  number  of 
these  ideas  have  been  applied  to  geophysical  problems  and  argue  for  their  continued 
refinement  and  application.  This  list  of  examples  is  representative  -ather  than  exhaus¬ 
tive.  In  particular  we  focus  on  the  attractor  reconstruction  technique  because  it  provides 
a  way  to  extract  the  greatest  amount  of  modeling  information  from  observational  data. 
The  attractor  that  is  reconstructed  from  the  data  is  shown  to  clearly  distinguish  between 
noise  and  chaos,  and  since  the  way  to  influence  systems  contaminated  by  noise  are  quite 
different  from  those  manifesting  fluctuations  due  to  low  order  nonlinear  interactions, 
being  able  to  distinguish  between  the  two  is  often  crucial.  When  such  an  attractor  can  be 
reconstructed  from  a  time  series  it  explicitly  gives  the  number  of  variables  required  to 
faithfully  model  the  phenomenon  of  interest. 

In  this  section,  our  discussion  spans  the  realm  of  activity  from  that  of  climate  to 
weather  to  fully  developed  wave  fields  on  the  ocean  surface. 

3.1  Weather  and  Climate  Attractors 

A  dominant  characteristic  of  meteorological  data  is  its  extreme  variability.  As  dis¬ 
cussed  by  Monin  (1972),  it  is  this  broadband  response  that  makes  the  predictability  of 
weather  patterns  from  deterministic  primitive  equations  so  uncertain.  The  fluctuations  in 
atmospheric  flow  field  that  give  rise  to  indeterminacy  have  been  associated  with  small- 
i  scale  turbulence.  The  scales  of  these  fluctuations  are  unresolved  in  global  circulation 

models  even  though  their  effects  are  manifest  on  these  large  scales.  One  technique  for 
treating  these  fluctuations  analytically  is  to  replace  the  deterministic  equations  by  sto¬ 
chastic  equations  [eg.  Thompson,  1957;  Landau  and  Lifshitz,  1959;  Lorenz,  1969; 
1  Hasselman,  1976;  Egger,  1982].  Such  replacements  have  in  the  past  been  purely 

phenomenological,  but  more  recent  studies  [West,  1982;  Lindenberg  and  West,  1984] 
suggest  how  one  may  proceed  more  systematically  from  the  deterministic  to  the  stochas- 
(  tic  representation.  A  third  alternative  has  been  to  replace  these  complex  field  models 

with  low-order  dynamical  models  having  chaotic  solutions,  such  as  first  done  by  Lorenz 
(1963).  It  is  this  last  approach  that  we  adopt  in  this  section. 
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We  discussed  the  Lorenz  model  at  some  length  in  Section  (2.1b)  and  showed  that 
the  basic  structure  of  the  observed  flow  depends  on  only  a  small  number  of  degrees  of 
freedom  in  phase  space.  The  general  hypothesis  (hope)  is  that  this  may  be  true  for 
observed  flow  fields.  Nicolis  and  Nicolis  (1984)  were  the  first  to  apply  the  attractor 
reconstruction  technique  to  a  geophysical  data  set  in  hopes  of  finding  a  climatic  attractor. 
They  applied  the  technique  to  the  isotope  record  of  a  deep  sea  core  to  estimate  the 
number  of  variables  governing  the  long-time  climate  evolution  during  the  past  million 
years.  They  used  the  method  of  Grassberger  and  Procaccia  (1983)  to  measure  the  dimen¬ 
sionality  of  the  data  set,  and  thereby  deduce  the  number  of  degrees  of  freedom  necessary 
to  control  the  underlying  dynamics  of  the  climate.  They  used  500  single  variable  values 
of  the  isotope  record  which  were  actually  interpolated  from  184  actual  measurements. 
The  value  of  the  correlation  dimension  they  obtained,  v  =  3.1,  was  criticized  by 
Grassberger  (1936)  as  resulting  from  the  smoothing  process  implicit  in  extending  the 
data  set  and  not  from  the  dynamics  of  the  climate  attractor.  Subsequent  analyses  have 
not  suffered  from  such  criticism  [Essex,  Lookman  and  Nerenberg,  1987;  Tsonis  and  Eis¬ 
ner,  1988;  Fraedrich,  1986]. 

In  Section  (2.3a)  we  discussed  the  correlation  dimension  as  a  measure  of  the  kind  of 
attractor  on  which  the  dynamics  for  a  given  data  set  exists.  Here  we  follow  the  more 
intuitive  approach  used  by  Fraedrich  (1986)  toward  a  general  definition  of  dimension.  A 
number  of  small  boxes  N(L),  each  of  side  L ,  are  used  to  fill  a  d  -dimensional  volume  V : 

V  =  £  Ld  =N(L)Ld  .  (3.1.1) 

N(L) 


Thus,  for  a  constant  volume  V,  the  number  of  boxes  filling  the  volume  increases 
inversely  with  L~d .  The  dimension  is  obtained  from  (3.1.1)  to  be 


In  fii(L)  In  V 

InL  +  InL  ’ 


(3.1.2) 


Taking  the  limit  of  decreasing  box  size  leads  to  the  more  general  definition  of  dimension 


i.taMil 

L  -*o  In  ( 1/L ) 


(3.1.3) 


first  developed  by  Hausdorff  and  more  recently  called  the  fractal  dimension  by  Mandel¬ 
brot  (1977).  If  the  volume  V  denotes  the  attractor  in  phase  space,  then  N(L)  is  the 
number  of  boxes  required  to  cover  the  attractor,  and  the  attractor  has  a  fractal  dimension 
d. 
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The  exponent  d  determines  how  fast  N (L)  increases  with  decreasing  L  and  may  be 
determined  by  plotting  in A/(L)  versus  ln(l/L).  However,  algorithms  to  count  these 
boxes  have  been  shown  to  be  impractical  in  determining  the  dimension  of  chaotic  attrac¬ 
tors  because  of  their  slow  convergence.  This  motivated  Grassberger  and  Procaccia 
(1983)  to  develop  another  method  which  would  estimate  a  lower  bound  on  the  fractal 
dimension,  i.e.,  the  correlation  dimension  in  Section  (2.3a).  Recall  that  X(  tj )  denotes  the 
trajectory  at  time  tj  in  an  m -dimensional  phase  space  X  =  (X(1),X(2),  •  ■  •  The 

number  N  (  )  of  pairs  of  points  whose  separation  distance  is  smaller  than  r , 
r  >|  X(  tj )  -  X(  tk )|  ,  is  formally  determined  by 

N(r)=  £  0[r-|X(r,)-X(f*)|]  (3  14) 

},k  =  \ 

and  N  is  the  total  number  of  points.  The  cumulative  distribution  function  C(r)  is  nor¬ 
malized  by  the  total  of  N2[~N(N  -  1)]  pairs  and  describes  how  the  number  of  pairs 
grows  with  increasing  distance  threshold  r .  For  N  — »«>,  the  growth  rate  changing  with 
the  dimension  d  is  determined  by 

C(r)=  lim^p-r*  .  (3.1.5) 

For  example,  consider  data  points  homogeneously  distributed  on  a  line,  then  the  number 
of  all  points  that  are  up  to  a  distance  r  apart  grow  linearly  with  r ,  i.e.,  are  proportional  to 
r .  For  data  homogeneously  distributed  over  a  plane  (in  a  volume),  the  number  of  points 
grows  quadratically  (cubicaliy)  with  r,  i.e.,  are  proportional  to  r2(r3).  Thus,  the  dimen¬ 
sion  of  the  attractor  in  phase  space  can  be  obtained  from 

d=\n C(r)_  (3J6) 

Inr 


where  C(r)  is  the  frequency  distribution  of  distance  of  the  pairs  of  points  situated  on  the 
time  trajectory  of  the  dynamical  system  in  an  m  -dimensional  phase  space. 

Because  we  do  not  know  the  dimension  of  the  attractor  in  advance,  one  successively 
determines  the  correlation  dimension  d  and  correlation  integral  C(r),  for  various  values 
of  the  embedding  dimension  m ,  i.e.,  d ( m )  and  Cm  ( r ): 


d(m) 


lnCm(r) 
In  r 


(3.1.7) 
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and  d(m)-*v  as  m  increases.  This  corresponds  to  the  m  -dimensional  phase  space 
X(t)=  {X(t),X(t  +  z),  ■  ■  ■  ,X(t  -  [m  -  1  ]x)}  where  T  is  the  time  lag  [Ruelle,  1981]. 
In  experimental  data,  noise  always  destroys  part  of  the  attractor.  If  the  noise  has  a 
characteristic  scale  rnoise ,  then  in  an  m  -dimensional  phase  space  the  noisy  trajectory  is 
space  filling  on  a  length  scale  r  smaller  than  rnoise : 

cm(r)~rm  for  r<rnoise  .  (3.1.8) 

Whereas  for  length  scales  greater  than  that  of  the  noise  level: 

Cm(r)~rd  for  r  ^  r noise  *  (3-19) 

Thus,  on  a  plot  of  In C(r)  versus  lnr,  there  should  be  a  break  in  the  curve  denoting  the 
change  in  slope  between  (3.1.8)  and  (3.1.9).  The  position  of  the  break  r  =  rMise  supplies 
information  on  the  noise  level  of  the  system. 

(a)  Surface  Pressure  and  Relative  Sunshine 

To  test  the  above  ideas  on  real  data  sets  Fraedrich  (1986)  selected  daily  values  of 
surface  pressure  and  relative  sunshine  duration  as  time  series.  In  Figure  (3.1.1)  is  dep¬ 
icted  a  three  dimensional  embedding  of  the  single  pressure  time  series  p  ( r )  for  three  dif¬ 
ferent  time  delay.  Note  how  the  trajectory  changes  from  one  tightly  constrained  to  nearly 
a  plane  when  T  =  3  hours,  to  one  that  seems  to  fill  the  three  dimensional  volume  for  z  =  3 
days. 

In  Figure  (3.1.2)  the  In  Cm(r)  versus  lnr  is  given  for  different  embedding  dimen¬ 
sions  ( 1  <m  <20)  for  z  -  3  days  corresponding  to  the  decorrelation  time  of  the  synoptic 
disturbances.  The  cumulative  distribution  of  pressure  distances  grows  with  increasing 
threshold  distance  r .  When  r  reaches  its  upper  limit,  the  distribution  function  converges 
to  unity.  Distribution  functions  are  shown  for  15  year  records  and  the  related  random 
series  having  the  same  mean,  variance  and  number  of  data  points.  A  random  data  set  of 
finite  length  (number  of  observations),  and  produced  by  a  random  number  generator,  is 
not  expected  to  follow  the  proportionality  d  =  m,  but  merely  the  inequality  d<m.  In 
Figure  (3.1.2)  we  see  that  the  observational  data  and  the  random  time  series  for  the  daily 
surface  pressure  appear  quite  similar.  In  Figure  (3.1.3)  Fraedrich  shows  the  d{m )  versus 
m  for  a  number  of  time  lags  and  we  see  that  there  is  some  tendency  for  d  ( m )  to  saturate 
at  a  value  d„,  but  what  that  value  is,  is  not  clear  from  the  figure.  It  is  clear,  however,  that 
the  dimension  is  below  that  of  the  random  time  series.  Fraedrich  points  out  that  the  data 


* 


•J 


495.chapter.3 


9-13-'88 


in  Figure  (3.1.3)  may  not  be  sufficiently  embedded  uecause  the  time  series  include 
weather  phenomena  from  winter  and  summer,  the  long-range  process  from  season  to  sea¬ 
son  and  the  interannual  variability. 

Separating  the  records  into  winter  seasons  (commencing  on  November  1)  and  sum¬ 
mer  seasons  (commencing  on  May  1),  Fraedrich  processes  the  data  for  14  winters  and  15 
summers  in  Figures  (3.1.3b)  and  (3.1.3c).  In  these  figures  the  dimension  values  clearly 
saturate;  d^- 3.2  for  the  winter  seasons  and  d^- 3.9  for  the  summer  seasons.  Thus  we 
conclude  that  although  summers  and  winters  separately  have  a  weather  attractor,  that  the 
combined  data  does  not.  Such  a  situation  could  arise  if  there  are  distinct  flow  regimes  in 
phase  space  for  each  of  the  seasons  rather  than  just  one  region  for  which  the  value  of  a 
parameter  is  changing  between  seasons. 

Similar  results  were  obtained  by  Fraedrich  for  the  sunshine  duration  data,  d^~ 3.1 
for  winter  and  dM~ 4.3  for  summer,  and  the  zonal  wave  amplitudes  of  the  500  mb  geopo¬ 
tential  waves  along  50  °  N  (dm~ 3.0  for  winter  and  d^~ 3.6  for  summer). 

Essex,  Lookman  and  Nerenberg  (1987),  as  well  as  Fraedrich,  also  used  these  ideas 
from  dynamical  systems  theory  for  the  study  of  global  climate.  These  authors  used  two 
variants  of  the  embedding  method  on  nine  sets  of  12,084  measurements  yielding  108,756 
values  in  all  of  daily  local  geopotential  values  at  500  mb  taken  at  12  UT  extending  over 
a  span  of  40  years.  These  data  are  very  like  one  of  the  sets  used  by  Fraedrich,  however, 
the  earlier  study  had  less  than  7.1  x  103  data  points.  The  first  method  of  Essex  et  al., 
treated  the  data  from  each  location  as  independent  measurements  from  a  single  site  by 
concatenating  the  different  time  series  from  each  site  to  create  one  large  time  series.  The 
results  of  this  method  are  shown  in  Figure  (3.1.4)  and  summarized  by  the  crosses  in  Fig¬ 
ure  (3. 1 .5).  The  straight  line  fits  in  Figure  (3. 1 .4)  of  In  Cm  ( r )  to  In  r  are  fairly  good  over 
a  reasonable  range  of  r . 

The  second  method  used  the  data  from  each  site  as  a  separate  coordinate  of  a  point. 
The  embedding  was  done  by  introducing  the  data  from  progressively  more  sites  as  addi¬ 
tional  coordinates.  The  results  are  indicated  in  Figure  (3.1.5)  as  triangles.  Figure  (3.1.6) 
shows  the  value  of  d(m )  at  each  r  for  m  =  9  using  this  method.  In  spite  of  the  variation 
in  the  data  points  the  plateau  region  is  quite  clear.  In  Figure  (3.1.5)  both  methods  are 
seen  to  yield  a  convergence  to  d„~ 6.  This  suggest  that  on  the  time  scale  of  decades,  cli¬ 
mates  might  be  represented  as  a  system  with  as  few  as  seven  degrees  of  freedom.  Thus 
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Essex  et  al.,  agree  with  the  earlier  conclusion  of  Nicolis  and  Nicolis  that  the  atmosphere 
and  oceans  exhibit  properties  of  a  low  dimensional  attractor,  although  the  two  analysis 
involve  data  that  differ  by  orders  of  magnitude. 

This  shortening  of  the  time  scales  was  brought  to  its  logical  conclusion  by  Tsonis 
and  Eisner  (1988)  who  used  data  on  the  order  of  hours,  again  orders  of  magnitude  shorter 
time  scales  than  those  considered  by  Essex  et  al..  Nicolis  and  Nicolis  using  single¬ 
variable  values  of  the  oxygen  isotope  records  of  deep-sea  cores  spanning  the  past  million 
years  calculated  a  climatic  attractors  dimension  of  3.1.  Fraedrich  and  Essex  et  al., 
analysed  weather  data  over  time  scales  of  15  to  40  years  and  obtained  a  dimension  of  the 
weather  attractor  of  between  6  and  7.  Tsonis  and  Eisner  using  wind  velocity  measure¬ 
ments  over  time  scales  of  minutes  calculated  a  weather  attractor  of  7.3. 

The  time  series  consisting  of  10  second  averages  of  the  vertical  wind  velocity 
recorded  10  m  above  the  ground  over  an  1 1  hour  period  is  shown  in  Figure  3.1.7a.  The 
total  number  of  data  points  is  3,960.  The  autocorrelation  function  and  spectral  density 
are  also  in  Figure  (3.1.7).  From  the  autocorrelation  function  a  conservative  estimate  of 
the  time  lag  was  made;  20  seconds.  In  Figure  (3.1.8)  the  logarithm  of  the  number  of 
pairs  versus  In  r  for  x  =  10 sec  for  embedding  dimensions  m  -  4, 6, 8, 10, 12  are  depicted. 
The  scaling  regions  are  indicated  by  the  straight  line  segments  indicated.  These  d(m)  are 
denoted  by  crosses  in  Figure  (3.1.9)  where  as  m  increases  we  see  that  d(m)—>d00  =  7.3. 
Therefore,  we  can  conclude  that  the  system  reresented  by  the  vertical  wind  velocity 
series  possesses  an  attractor.  The  non-integer  dimension  of  the  attractor  suggests  that  it 
is  chaotic,  i.e.,  that  the  submanifold  is  a  fractal  set. 

3.2  Fractal  Dimension  and  the  Ocean  Surface 

Recently  the  importance  of  nondifferentiable  curves  and  surfaces  have  become 
apparent  in  different  areas  of  the  physical  science.  This  growing  awareness  is  due  in  no 
small  part  to  the  efforts  of  Mandelbrot  (1977,  1982)  in  drawing  attention  to  the  fractal 
paradigm  in  modeling  a  wide  variety  of  natural  phenomena.  We  have  pointed  out  that  a 
fractal  object  is  one  that  possesses  no  smallest  scale.  Of  course  in  real  systems  there  is 
always  a  largest  and  smallest  scale,  it  is  a  question  of  the  size  of  the  interval  overwhich 
no  characteristic  scale  is  apparent.  For  some  applications  a  single  decade  of  data  is 
sufficient,  for  others  one  would  want  two  or  three  decades  of  scale-free  data  before  the 
process  was  pronounced  a  fractal.  In  a  geophysical  context  Ausloos  and  Berman  (1985) 
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have  recently  used  a  fractal  function  to  model  undersea  topography  and  Stiassne  (1988) 
used  one  to  model  the  surface  of  the  ocean.  We  review  these  two  applications  in  this  sec¬ 
tion. 

The  German  mathematician  Karl  Weirerstrass  cast  the  argument  for  fractals  into  a 
particular  mathematical  form.  His  intent  was  to  construct  a  series  representation  of  a 
continuous,  non-differentiable  function.  His  function  was  a  superposition  of  harmonic 
terms: 

W(t)=  2  — cos(6na>0O  (3.2.1) 

n=  oa" 

a  fundamental  with  frequency  Wq  and  unit  amplitude,  a  second  periodic  term  of  fre¬ 
quency  b  co0  with  amplitude  1  la ,  a  third  periodic  term  of  frequency  b2 co0  with  amplitude 
1/a2,  and  so  on.  The  resulting  function  is  an  infinite  series  of  periodic  terms,  each  term 
of  which  has  a  frequency  that  is  a  factor  b  larger  (&>1)  than  the  preceding  term  and  an 
amplitude  that  is  a  factor  1/a  smaller  (a  >1). 

Note  that  for  this  concept  of  a  fractal  function,  or  fractal  set,  there  is  no  smallest  or 
characteristic  scale.  For  b  >  1  in  the  limit  of  n  terms,  n  — the  frequency  bn co0— »°°,  and 
there  is  no  highest  frequency  contribution  to  the  Weierstrass  function.  Of  course,  if  one 
thinks  in  terms  of  periods  rather  than  frequencies,  then  the  shortest  period  contributing  to 
the  series  is  zero.  Consider  what  is  implied  by  this  lack  of  a  smallest  scale  in  period,  or 
equivalently  by  the  lack  of  a  largest  scale  in  frequency.  Imagine  a  continuous  line  on  a 
two-dimensional  Euclidean  plane  and  suppose  that  the  line  has  a  fractal  dimension 
greater  than  unity  but  less  than  two.  How  would  such  a  curve  appear?  In  this  case  we 
are  superimposing  smaller  and  smaller  wiggles:  the  curve  would  be  like  the  irregular  line 
on  a  map  representing  a  very  rugged  sea  coast. 

At  first  glance  this  curve  would  seem  to  be  a  ragged  line  with  many  abrupt  changes 
in  direction,  as  in  Figure  (3.2.1a).  If  we  now  magnify  a  small  region  of  the  line,  as  in 
Figure  (3.2.1b),  we  see  that  the  enlarged  region  appears  qualitatively  the  same  as  the  ori¬ 
ginal  curve.  If  we  now  magnify  a  small  region  of  this  new  line,  as  in  Figure  (3.2.1c),  we 
again  obtain  a  curve  that  is  qualitatively  indistinguishable  from  the  first  two.  The  pro¬ 
cedure  can  be  repeated  indefinitely  [as  long  as  n  -*<»  in  (3.2.1)],  and  the  behavior  of  the 
mathematical  function  is  analogous  to  the  discontinuous,  inhomogeneous  clumping  of 
matter  in  space.  The  equivalence  of  the  function  from  one  scale  to  the  next  is  called 
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self-similarity  and  the  measure  of  the  degree  of  self-similarity  in  the  Weierstrass  function 
(in  terms  of  frequency  and  amplitude),  is  precisely  the  fractal,  or  Hausdorff-Bessiovich 
dimension:  d  =\na/\nb. 

Because  of  the  infinite  layers  of  detail,  one  cannot  draw  a  tangent  to  a  fractal  curve, 
which  means  that  the  function,  although  continuous,  is  not  differentiable.  It  also  means 
that  the  more  terms  one  allows  to  contribute  to  the  function  (the  more  accurately  one 
maps  the  coastline),  the  longer  the  curve  becomes,  and  there  is  no  upper  limit  to  its 
length. 

Ldvy  introduced  a  generalization  of  the  Weierstrass  function  which  was  later  dis¬ 
cussed  by  Mandelbrot  (1977)  and  subsequently  analysed  in  extensive  detail  by  Berry  and 
Lewis  (1980).  This  function  was  called  the  Weierstrass-Mandelbrot  function  by  the  last 
authors,  but  here  we  refer  to  it  as  the  extended  Weierstrass  function  since  there  will  be  a 
number  of  such  generalizations.  Here  again  we  have  a  superposition  of  sinusoidal  terms 
with  geometrically  spaced  frequencies  and  amplitudes  that  follow  an  inverse  power  law: 

W(t)  =  £  b-n{2-D)[\-eib'Wot]ei<3'  ,  (3.2.2) 

n  -  -  oo 

where  the  phase  <(>„  is  arbitrary,  \<D  <  2,  and  we  have  set  2  -  D  =  In  a  /In  b .  In  the  fol¬ 
lowing  discussion  we  scale  time  in  such  a  way  that  w0  =  1 .  For  a  complete  discussion  of 
the  properties  of  (3.2.2)  we  refer  the  reader  to  Berry  and  Lewis  (1980);  here  we  merely 

►  record  some  of  them  that  are  important  for  our  discussion. 

The  set  of  phases  {<J>„  }  may  be  chosen  either  deterministically  or  randomly.  If  <}>„  is 
a  random  variable  uniformly  distributed  on  the  interval  (0,  2k),  then  each  choice  of  {<>„  } 
constitutes  a  member  of  an  ensemble  for  the  stochastic  function  W(t).  If  the  phases  are 

t 

also  independent  and  b— >1  + ,  then  W(r)  is  a  Gaussian  random  function.  The  condition 
\<D  <2  is  required  to  ensure  the  convergence  of  the  sum. 

To  generalize  the  extended  Weierstrass  function  Ausloos  and  Berman  preserved  two 
i  properties:  homogeneity  and  scaling  (statistical  properties).  Consider  the  increments  of 

W: 

&W(t,x)  =  W(t  +x)- W(f) 

>  =  £  b~n{2~D)[eib't ,  (3.2.3) 

«=—<*> 

and  assume  that  the  4>n  are  independent,  random  variables  uniformly  distributed  on  the 


» 
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interval  (0,  2k).  The  mean  square  increment  is  then 
V(x)  =  <|  AW(r,x)|  2> 

=  £  b-n(4~2D)2[\-cos(bnx)]  (3.2.4) 

n  =  -  <*» 


where  the  brackets  denote  an  average  over  an  ensemble  of  realizations  of  the  <J)„- 
fluctuations.  The  right  hand  side  of  (3.2.4)  is  independent  of  t,  i.e.,  it  depends  only  on 
the  time  difference  x,  so  that  W  ( t )  is  homogeneous  (also  called  stationary  when  t  refers 
to  the  time). 

Equation  (3.2.4)  also  leads  to  the  scaling  property  of  the  extended  Weierstrass  func¬ 
tion: 

V(fcx)  =  h2(2_D)V(x)  .  (3.2.5) 


If  the  sum  giving  the  extended  Weierstrass  function  (3.2.2)  covered  the  interval  from 
n  =  0  to  n  =  »,  then  the  scaling  property  expressed  in  (3.2.5)  would  hold  only  approxi¬ 
mately.  A  two-variable  scalar  extended  Weirerstrass  function  W(  r)  is  proposed  by 
Ausloos  and  Berman  to  be  one  such  that 

V(  p)  =  < |  AW ( r,  p)  |  2  >  (3.2.6) 


so  that  for  all  r 

V(bp)  =  b2(3~D)V(p )  (3.2.7) 

where  in  this  case  2<D  <  3. 


After  discussing  some  candidate  functions  that  for  one  physical  reason  or  another 
were  considered  unsatisfactory,  Ausloos  and  Berman  decided  on  one  that  is  a  random 
superposition  of  weighted  “ridge-like”  surfaces: 

</2 


W(r)  = 


In  b 


M 


L  Arn  Z  (k0bn)D-i[l-e‘b'korcos(Q  ot")] e<6m"  .(3.2.8) 

m  =  1  n  =  —  oo 


Here  is  a  wave  number  that  can  be  used  to  scale  horizontal  variations;  the  normaliza¬ 
tion  factor  is  chosen  to  make  the  series  converge  as  b— >\  +  and  M  — the  angle  am 
gives  the  orientation  of  the  corrugation  of  the  surface  having  an  amplitude  Am ;  and  the 
phases  are  again  defined  as  random  variables.  It  is  a  simple  matter  to  verify  that 
(3.2.8)  has  the  scaling  property  of  (3.2.7)  and  that  V(  p)  is  homogeneous. 
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Ausloos  and  Berman  present  a  number  of  computer  plots  of  the  surfaces  generated 
by  the  extended  Weirerstrass  function  (3.2.8).  Their  interest  was  in  measuring  the  effect 
of  the  small  length  scales  on  long-range  acoustic  propagation  in  the  ocean.  Here  we  are 
more  interested  in  the  possible  use  of  this  function  to  model  the  ocean  surface.  Therefore 
let  us  examine  some  of  their  computer  plots  interpreted  as  snap-shots  of  the  sea  surface. 

In  Figure  (3.2.2)  we  see  a  figure  that  could  be  interpreted  as  a  one-dimensional  sea 
surface  consisting  of  long  ridges  of  random  height  and  wavelength.  This  is  not  ’nlike  the 
surface  profiles  shown  on  pages  330  and  331  of  Neumann  and  Pierson  (1966).  In  Figure 
(3.2.3)  two  of  these  surfaces  at  right  angles  are  superposed  to  give  rise  to  a  two- 
dimensional  random  surface  which  resembles  that  of  a  high  sea.  To  the  unaided  eye  the 
most  convincing  example  of  a  sea  surface  is  given  in  Figure  (3.2.4).  Here  we  seem  to  see 
the  ever  present  small  scale  structure  riding  on  top  of  a  large  scale  wave. 

A  different  extended  Weierstrass  function,  but  one  that  has  both  properties  of 
(3.2.1)  and  its  extension  to  two  spatial  dimensions  (3.2.2)  is  considered  by  Falconer 
(1985) to  be 

fl  nh}*  M 

W(r)=  — -  2  Am  2)  b  n(3  D)cos[*oi>n(xcos9m  +  jysin9m)-f  .  (3.2.9) 

„  M  )  m  =  1  n  =0 

If  we  interpret  the  phase  as  the  time-dependent  quantity 


§mn  ^ 


(3.2.10) 


where  co„  is  the  frequency  of  a  deep  water  gravity  wave  having  wavenumber  kn  =  k0bn , 

^n^<Wn=<gQ>n/2  (3-2.11) 

and  g  is  the  acceleration  of  gravity.  Thus  we  have 

Wr(r)=  2  Am  £  i-^-^cos  [^"(jccos^  +ysin0m) 

M  m  =  1  n  =0 


+  4 >mn  ~0>o bn,2t)  , 


(3.2.12) 


in  which  a  superposition  of  propagating  waves  has  been  made  to  form  W  ( r).  Following 
Stiassnie,  the  amplitudes  Am  can  be  chosen  in  such  a  way  that  the  mean  square  value  of 
W{  r)  is  the  same  as  that  of  the  sea  surface.  To  do  this  we  recall  that  the  surface  wave 
spectrum  is  given  by 


Bu*  n 

vP(k)=  —  G(9)/*p 

.  8 
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where  the  function  G(0)  accounts  for  the  directionality  of  the  waves,  u,  is  the  wind 
function  velocity  and  B  =0.012.  Noting  that  we  have  selected  a  geometric  series  in 
wavenumber  we  replace  the  differential  kdkdQ,  by  the  discrete  versions  0m  =  mA0  and 
dkn  =kn\nb  so  that 

¥( k)kdkd 0-3-*„ dkn A2  =  — — ^-ln bG{Qm)k2~^2~^  .  (3.2.14) 


Using  (3.2.14)  to  replace  Am  in  (3.2.12)  therefore  yields 


2w  Bu *  In  b 

W(r)  =  — - 

\  *  f  >  m  I/.,  'l/' i  r 


M  gHr-D)\  nX  1 


IYI  OO 

S  [G(0m)]*  I  i»-n(3-£>) 


cos  [£o£n(rcos0m  +  y sin0m )  +  6m*  -(0Qbn/2t]  (3.2.15) 


if  the  exponents  of  k0  are  equal,  i.e., 

D  -  4-  (3/2  (3.2.16) 

Thus  the  surface  represented  by  (3.2.15)  has  a  fractal  dimension  D  dependent  on  the 
index  of  the  observed  surface  wave  spectrum.  If  we  take  (3  =  7/2  as  suggested  by 
Kitaigorodskii  (1983)  then  the  fractal  dimension  of  the  sea  surface  is  D  -  2.25. 

As  Stiassne  points  out,  surface  tension  introduces  a  cut-off  at  say  k  =  ks  so  that  the 
actual  sea  surface  is  only  approximately  fractal.  It  has  a  fractal-like  structure,  over 
length  scales  from  2n/ks  to  2 K/k0,  which  corresponds  to  approximately  0.1  m  to  100  m  in 
the  open  ocean. 
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4.  Water  Wave  Attractor 

We  now  turn  to  the  question  of  whether  these  dynamic  techniques  are  useful  in 
interpreting  data  sets  from  either  the  open  ocean  or  from  wave  tanks.  Traditionally  we 
think  of  surface  wave  interaction  phenomena  as  consisting  of  both  broad  band  and  nar¬ 
row  band  processes.  The  mode/mode  energy  transfer  resulting  from  the  nonlinear  hydro- 
dynamic  interactions  is  broad  band,  as  is  the  general  wave  train  instability.  On  the  other 
hand,  the  Benjamin  Feir  instability  mechanism  and  recurrence,  as  well  as  surface 
envelope  solitons  are  narrow  band  phenomena.  These  processess  have  been  understood 
on  the  basis  of  the  weak  interaction  theory  of  interacting  nonlinear  gravity  waves  [see  eg. 

West  (1981)]. 

Many  of  these  weak  interaction  effects  are  discussed  by  West  (1981)  and  as  sum¬ 
marized  by  Phillips  (1Q?Q)  these  phenomena  share  the  following  characteristics: 

(1)  The  weak  interaction  phenomena  may  be  manifest  only  by  the  longer  (dom¬ 
inant)  waves  in  the  wave  field ;  the  behavior  of  short  waves  (wavelengths  less  than  the 
dominant  one)  is  determined  by  strong  interactions  with  long  waves  rather  than  by  weak 
interaction  processes.  This  is  due  in  part  to  the  fact  that  only  the  dominant  waves  main¬ 
tain  their  integrity  sufficiently  long  for  the  mechanism  to  work. 

(2)  They  are  evolutionary  phenomena  with  a  common  time  scale  of  order  e~2  times 
the  wave  period,  where  e  is  the  slope  of  the  dominant  wave.  Under  most  oceanic  condi¬ 
tions,  this  time  scale  is  of  order  100  wave  periods  (minimally)  of  the  dominant  wave. 

Equivalently,  the  distance  that  waves  travel  before  these  effects  become  significant  is  of 
order  e-2  times  the  dominant  wavelength.  The  shorter  waves  disappear  and  are  refreshed 
by  the  wind  on  a  much  shorter  space-time  scale.  These  interactions  lead  to  Gaussian 
statistics  of  the  surface  wave  field  [Benny  and  Newell  (1967),  Hasselmann  (1968)]. 

(3)  The  state  of  a  wave  field  at  a  given  point  is,  in  part,  the  result  of  these  processes 

acting  over  the  preceding  time,  but  the  modification  they  can  produce  in  the  wave  field  in 
a  limited  distance  or  time  (a  few  wavelengths  or  periods)  is  small,  of  order  e-2,  character¬ 
istically  one  percent  or  less.  Consequently,  they  are  intrinsically  unable  account  for  any 
observed  local  response  in  the  wave  field  to  a  local  disturbance  such  as  the  modulation  of 
a  short  wave  near  the  peak  of  a  long  wave  in  a  wind  generated  field.  • 
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The  strong  interactions  are  those  which  occur  on  a  time  scale  of  the  same  order  as 
the  wave  period  of  interest.  They  include  the  instability  leading  to  wave  breaking,  wave 
breaking  itself,  wave-current  interactions,  the  distortion  of  a  short  wave  energy  packet  by 
a  long  wave,  the  process  of  parasitic  capillary  formation  and  micro-scale  breaking 
induced  by  surface  wind  drift.  These  phenomena  provide  the  possibility  of  a  rapid 
response  to  oceanic  perturbations  and  are  much  less  studied  than  the  above  weak  interac¬ 
tions.  Again  summarizing  the  discussion  of  Phillips  (1979)  and  West  (1981): 

(1)  The  surface  wave  train  modulation  produced  by  internal  waves  or  other  currents 
is  a  strong  interaction  effect  in  that  it  can  be  produced  by  a  single  internal  wave  pulse 
(though  it  is  frequently  described  as  a  “resonance,”  a  weak  interaction  phenomena). 

(2)  A  wave  train  or  packet  of  short  waves,  interacting  with  a  longer  wave,  varies  in 
amplitude  and  scale  with  respect  to  its  position  on  the  long  wave.  The  modulation  pat¬ 
tern  of  the  packet  has  a  spectral  signature  that  is  distributed  over  a  range  of  frequencies 
proportional  to  the  long  wave  speed  and  is  approximately  proportional  to  the  long  wave 
slope.  As  the  long  wave  slope  increases,  the  spectral  signature  is  dominated  more  and 
more  by  the  properties  of  the  short  waves  near  the  long  wave  crest ,  and  less  by  those  at 
other  points  of  the  long  wave  cycle.  This  implies  that  the  surface  fluctuations  cannot  be 
homogeneous  in  space. 

(3)  Micro-scale  surface  properties  such  as  the  formation  of  parasitic  capillary 
waves  and  microscale  breaking,  i.e.,  the  formation  of  sheets  of  bubbles,  have  densities 
(number  of  occurrences  per  unit  area)  that,  in  an  active  wind-generated  sea,  respond 
rapidly  to  energy  exchanges  with  short  gravity  waves  and  more  sensitively  than  changes 
in  the  short  wave  amplitudes  themselves.  In  this  way,  they  can  serve  as  indicators  to  per¬ 
turbations  on  the  scale  of  internal  waves  and  other  internal  flow-fields.  This  appears  as 
large  scale  and  long  time  correlations  in  the  surface  statistics. 

Here  we  wish  to  determine  if  the  surface  wave  field  can  be  described  by  a  low¬ 
dimensional  attractor,  which  in  the  present  context  would  be  viewed  as  a  strong  interac¬ 
tion.  To  do  this  we  need  to  analyze  a  “time”  series  characteristic  of  the  sea  surface.  For 
this  we  use  a  numerically  generated  surface  wave  field  using  the  new  numerical  method 
for  surface  hydrodynamics  developed  by  West,  Brueckner,  Janda,  Milder  and  Milton 
(1987). 
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4.1  Numerically  Generated  Data 


We  use  the  computer  code  of  West  et  al.,  (1987)  to  numerically  generate  a  one¬ 
dimensional  surface  displacement  C,(x,t).  For  this  purpose  we  use  512  wavenumber 
modes  in  units  such  that  k0  =  2k/L  =  1  and  the  acceleration  of  gravity  is  unity: 


-  Im 


OIjL  ..  r—  . 

n  -  1 


(4.1.1) 


The  initial  mode  amplitudes  (<3n(0)}  are  selected  so  as  to  have  a  Phillips  spectrum  k~ 3 
and  independent  random  phases  uniformly  distributed  on  the  interval  ( 0,  2ji).  The  equa¬ 
tions  of  motion  given  in  West  et  al.  are  then  integrated  forward  in  time  to  their  values 
{an{t)}  on  a  Cray  computer.  The  surface  displacement  is  then  used  in  two  ways.  First, 
at  a  fixed  time  t ,  the  frozen  surface  on  the  spatial  interval  ( 0,  2k)  is  segmented  into  M 
points  which  constitute  the  “time”  series.  Secondly,  at  a  fixed  spatial  point,  x ,  the  varia¬ 
tion  in  the  surface  elevation  in  time  is  treated  as  the  time  series.  The  latter  use  of  the  data 
is  what  one  would  record  with  a  wave  staff  located  at  x  measuring  the  passage  of  a  wave 
field. 


We  compare  the  above  data  with  a  completely  random  data  set  generated  by  the 
function 

/?(0=  2  cos  (^r +<!>„)  (4.1.2) 

n  -  1 


where  <!>„  is  a  random  variable  uniform  on  the  interval  ( 0, 27t)  and  N  is  the  number  of 
modes  in  the  sum.  The  variable  R  ( / )  becomes  a  zero-centered  Gaussian  random  vari¬ 
able  for  large  N  which  we  also  segment  into  M  points  for  our  random  time  series. 

The  “experimental”  data  in  method  1  are  given  by  the  set  {^j=^(xj,  r)}  where  we 
divide  the  continuous  spatial  curve  into  1024  equally  spaced  points.  We  use  these  data  to 
construct  the  set  of  variables  _x,  ■  ■  •  _(m  _  jyJ  in  an  m  -dimensional  embed¬ 

ding  space  and  apply  the  Grassberger-Procaccia  correlation  algorithm  to  determine  the 
dimensionality  of  the  data.  In  Figure  (4. 1 . 1 )  we  depict  In  Cm  ( r )  versus  In  r  for  a  number 
of  embedding  dimensions.  In  Figure  (4.1.2)  we  show  the  correlation  dimension  for  the 
data  set  using  two  lag  times  x  as  a  function  of  embedding  dimension.  The  random  data, 
partitioned  in  the  same  way,  yields  the  result  d{m)-m  as  expected.  The  surface  data, 
however,  clearly  show  a  tendency  to  saturate  for  both  T  =  8  and  x  =  15.  This  implies  that 
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allhough  the  surface  was  generated  using  512  degrees  of  freedom  (in  the  spectral  sense) 
as  few  as  five  or  six  degrees  of  freedom  may  be  adequate  to  describe  its  spatial  behavior. 
Note  that  in  method  1  our  data  describe  points  in  space,  not  points  in  time.  Therefore  the 
interpretation  of  these  results  is  not  as  straightforward  as  the  earlier  applications. 

We  note  the  kink  in  the  In  Cm  ( r )  versus  In  r  curve  in  Figure  (4.1.1)  for  large  m .  At 
first  It  was  not  clear  how  one  might  interpret  this  result,  but  analyses  due  to  Theiler  in  his 
Ph.D.  thesis  (1987)  suggests  that  the  kink  may  be  due  to  correlation  in  the  data  set.  In  his 
analysis,  Theiler  considers  a  stochastic  data  set  that  is  Gaussian  with  a  prescribed  correla¬ 
tion  time  a.  The  correlation  time  is  observed  to  produce  a  kink  [cf.  Figure  4.1.3]  just  of 
the  form  seen  in  Figure  (4.1.1)  and  to  result  in  the  correlational  dimension  not  increasing 
linearly  with  the  embedding  dimension.  This  observation  suggests  that  the  longest 
wavelength  waves  on  the  sea  surface  modulate  the  short  wavelength  waves,  and  this 
modulation  is  picked  up  as  a  correlation  of  the  small  scale  structure  at  the  higher  embed¬ 
ding  dimensions.  Theiler  developed  a  technique  to  suppress  this  correlation  which  we 
propose  to  test  in  the  present  context.  The  method  involved  deleting  the  n-  nearest 
neighbors  in  the  calculation  of  the  correlation  function,  but  we  will  not  pursue  that  tech¬ 
nique  further  here. 

We  now  turn  to  the  second  method  of  collecting  data  in  whi^h  we  record  the  height 
of  the  sea  surface  as  a  function  of  time  at  a  fixed  spatial  location.  We  segment  the  spatial 
interval  ( 0,  2k)  into  16  equi-distant  points.  At  each  point  we  record  the  surface  elevation 
for  an  interval  of  time  four  units  long.  This  corresponds  to  two  periods  of  the  longest 
wave  present  on  the  sea  surface.  Each  of  the  16  time  series  is  discretized  into  1024 
points,  yielding  a  total  of  16,384  data  points  if  we  were  to  concatenate  all  the  data  as 
done,  for  example,  by  Essex  et  al.,  (1987)  in  their  global  climate  study. 

We  concatenate  the  data  from  eight  adjacent  locations  to  define  a  single  time  series 
consisting  of  8,192  data  points.  In  method  1  we  picked  the  value  of  t  by  inspection.  In 
the  present  method  we  were  guided  in  our  choice  by  calculating  the  time  at  which  the 
auto  correlation  function  first  goes  through  zero.  The  preferred  technique  would  have 
been  to  calculate  the  first  minimum  in  the  mutual  information,  but  calculating  the  auto¬ 
correlation  function  was  faster.  The  z  value  determined  by  the  autocorrelation  function  is 
50,  so  we  selected  a  value  substantially  above  that,  i.e.,  t  =  80.  In  Figure  (4.1.4)  we  see 
that  there  are  two  distinct  scaling  regimes  for  the  data  above  m  =  8.  The  largest  regime 
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has  v=2  and  is  fairly  independent  of  the  embedding  dimension.  The  second  regime  is 
plotted  in  Figure  (4.1.5)  as  a  function  of  embedding  dimension  and  is  seen  to  saturate 
around  v  =  d.^-6. 
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5.  Conclusions  and  Speculations 

The  physicist-mathematician  Mars  Kac  was  fond  of  pointing  out  the  difference 
between  a  demonstration  and  a  proof.  He  maintained  that  a  demonstration  would  con¬ 
vince  a  reasonable  man  of  the  truth  of  a  mathematical  statement,  whereas  it  required  a 
proof  to  convince  a  mathematician.  Herein  we  have  demonstrated  the  utility  of  certain 
generic  concepts  from  nonlinear  dynamical  systems  theory  in  geophysics  and  physical 
oceanography.  We  have  reviewed  how  the  phase  space  reconstruction  technique  can  be 
used  to  determine  the  minimum  number  of  variables  necessary  to  describe  certain  geo¬ 
physical  phenomena.  Since  we  are  accustomed  to  describing  geophysical  flows  by  fields 
which  are  spatially  continuous  and  therefore  have  an  infinite  number  of  degrees  of  free¬ 
dom,  this  notion  of  representing  such  flows  by  low-dimensional  attractors  is  admittedly  a 
surprising  one,  but  not  so  surprising  as  being  able  to  determine  the  dimensions  of  the 
attractor,  and  thereby  the  number  of  variables  needed  to  describe  the  system  directly 
from  a  single  time  series  of  sufficient  length. 

Let  us  recap  the  attractor  (phase  space)  reconstruction  technique.  We  assume  the 
system  of  interest;  the  isotope  concentration  from  a  deep  ocean  core,  the  duration  of 
sunshine  as  a  function  of  time  at  a  measuring  station,  the  speed  of  the  wind  at  a  given 
location  and  the  height  of  the  sea  surface  measured  by  a  wave  staff,  can  each  be 
described  by  m  -variables,  where  m  is  large  but  unknown.  At  any  instant  of  time  there  is 
a  pont  X(r)  =  {X(1Vr),X(2)(r),...,X(m)(r)}  in  an  m -dimensional  phase  space  that  com¬ 
pletely  characterizes  the  system.  The  point  moves  around  as  the  system  evolves.  In  the 
geophysical  cases  discussed  herein  the  motion  had  previously  appeared  to  be  purely  ran¬ 
dom  so  that  one  must  distinguish  between  a  system  confined  to  a  chaotic  attractor  from 
one  driven  by  noise.  For  many  reasons  geophysical  observations  are  often  restricted  to 
the  output  of  a  single  detector  which  selects  one  of  the  m  components  of  the  system  for 
monitoring.  In  general  the  scientist  does  not  know  the  dimension  of  the  phase  space 
since  the  important  dynamic  variables  are  usually  unknown  and  therefore  as  much  infor¬ 
mation  as  possible  must  be  extracted  from  the  single  time  series,  X(1)(r)  say.  For 
sufficiently  long  delay  times  t  one  use  the  embedding  theorem  to  construct  the  sequence 
of  displaced  time  series  {X(1)(r),X(1)(r  -  x), . .  .  ,X(I)[f  -  (m  -  1  )x] }.  In  this  Tri¬ 
dimensional  phase  space  one  then  determines  if  the  system  dynamics  is  confined  to  a 
low-dimensional  attractor. 


495/chapterJ 


9- 1 3-'88 


-63- 


This  procedure  was  followed  in  the  examples  listed  in  Section  3  and  the 
Grassberger-Proccacia  method  of  determining  the  correlation  dimension  was  used  to 
determine  if  chaotic  attractors  exist  in  geophysical  data  sets.  Nicolis  and  Nicolis  (1984) 
were  the  first  to  find  a  climate  attractor  albeit  with  an  insufficient  data  base  they  found 
v  =  3.1.  Essex  et  al.  (1987)  and  independently  Fraedrich  (1986)  found  climate  attrac¬ 
tors  but  with  a  dimension  more  than  twice  that  found  by  Nicolis  and  Nicolis.  Using  daily 
values  of  surface  pressure  and  sunshine  duration  Fraedrich  found  separate  weather  attrac¬ 
tors  for  summer  and  winter  with  dimensions  in  the  intervals  3.2  <  v  <  4.3.  Essex  et  al. 
used  the  local  geopotential  value  to  find  a  weather  attractor  with  v  =  6.  Tsonis  and  Eis¬ 
ner  used  vertical  wind  velocity  measurements  to  calculate  a  weather  attractor  with 
v  -  7.3.  The  time  series  used  to  determine  these  various  attractors  spanned  millions  of 
years  for  the  single  values  of  the  oxygen  isotope  records  of  deep-sea  cores  down  to  wind 
velocity  measurements  over  time  scales  of  a  few  hours. 

Our  own  calculation  of  water  wave  attractors  in  Section  4  suggests  that  the  one¬ 
dimensional  sea  surface  may  be  describable  by  as  few  as  two  degrees  of  freedom  in  one 
domain  and  as  many  as  six  in  another.  The  first  domain  is  not  unlike  the  result  obtained 
using  the  generalized  Weierstrass  function  to  model  the  sea  surface.  The  possible  rela¬ 
tion  between  these  two  approaches  has  not  as  yet  been  explored  in  the  present  geophysi¬ 
cal  context.  It  is  of  interest,  however,  to  consider  how  these  techniques  may  be  used  for 
the  remote  sensing  of  the  sea  surface. 

The  fundamental  problem  of  remotely  sensing  patterns  on  the  sea  surface  is  the 
detection  of  a  signal  in  a  large  dynamic  range  background.  This  can  be  seen  from  even  a 
superficial  examination  of  the  variety  of  physical  processes  underlying  the  signal  pro¬ 
cessing  problems  of  interest.  Consider  the  surface  of  the  ocean;  the  nonlinear,  wind¬ 
generated,  surface  disturbances  constitute  the  noisy  background,  [see  e.g.  West,  1981] 
and  the  modulation  of  them  by  surface  current  disturbances  constitute  the  signal.  Prob¬ 
ing  the  ocean’s  surface  via  radar  relies  on  an  understanding  of  the  ocean  surface  modula¬ 
tion  mechanism  since  the  scattering  is  from  very  short  Pragg  waves  (wavelengths  of  a 
few  centimeters)  and  the  currents  of  interest  have  scales  of  tens  of  hundreds  of  meters. 
Going  now  into  the  ocean  interior,  the  acoustic  signal  generated  by  rainfall  or  the  break¬ 
ing  of  surface  is  also  of  interest.  In  the  far  field  region  the  pressure  and/or  velocities  gen¬ 
erated  by  these  sources  have  been  contaminated  by  the  dynamics  of  the  intervening 
medium,  e.g.,  by  internal  waves,  patches  of  turbulence  and  acoustic  noise  from  the  sea 


495/chapter.5 


9-i*-'88 


-64- 


surface.  These  undersea  phenomena  also  influence  surface  currents  and  act  to  degrade 
the  signal.  The  motion  of  conducting  sea  water  also  perturbs  the  earth’s  magnetic  field, 
e.g.  internal  waves  give  rise  to  a  clear  magnetic  field  signature. 

The  remote  detection  and  identification  of  undersea  phenomena  is  therefore  a  prob¬ 
lem  of  indirect  measurement  in  that  one  must  be  able  to  correctly  interpret  the  distur¬ 
bances  of  the  ambient  surface  waves  using  radar  (or  light)  backscattered  from  the  sea  sur¬ 
face  and  the  perturbation  in  the  earth’s  magnetic  field  produced  by  the  currents’  passage. 

Techniques  for  processing  the  perturbed  radiation  whether  electromagnetic,  acoustic 
or  thermal  run  the  gamut  from  the  application  of  statistical  pattern  recognition  methods, 
to  procedures  for  coherent  data  acquisition,  to  artificial  neural  systems.  These  and  other 
techniques  have  their  proponents,  but  they  by-and-large  miss  an  essential  feature  of  the 
geophysical  problem;  the  sea  surface  and  interior  are  described  by  nonlinear  dynamic 
(hydrodynamics)  processes,  [Phillips,  1977], 

Historically  the  aperiodic  component  of  experimental  time  series,  such  as  the  back 
scattered  radiation,  has  been  interpreted  as  random  noise.  Thus  the  processing  pro¬ 
cedures  have  relied  on  techniques  that  separate  the  ‘signal’  from  the  ‘noise’  in  a  linear 
additive  way.  This  view  is  consistent  with,  and  is  probably  a  consequence  of,  the  tradi¬ 
tional  interpretation  of  stochastic  processes  in  the  physical  sciences.  It  results  from  a 
dynamic  system  being  coupled  to  a  complicated  environment,  i.e.,  the  environment 
induces  a  random  component  into  the  system  dynamics.  A  Fourier  decomposition  of 
such  a  time  series,  one  consisting  of  a  deterministic  signal  linearly  superimposed  on  a 
noise,  would  yield  a  spectrum  consisting  of  sharp  spectral  peaks  at  the  natural  frequen¬ 
cies  of  the  signal  and  a  broad  band  background  depicting  the  power  of  the  noise.  This 
interpretation  has  become  so  ingrained  that  its  inverse  is  often  applied,  that  is  to  say  that 
a  spectrum  such  as  shown  in  Figure  2.3.3,  consisting  of  a  broad  band  spectral  level  on 
which  a  number  of  peaks  are  superimposed,  often  suggests  the  linear  additive  model  of 
signal  and  noise.  That  interpretation  cannot  in  fact  be  applied  to  the  spectrum  in  Fig¬ 
ure  2.3.3  because  the  underlying  time  series  was  generated  by  a  nonlinear  equation  of 
motion  having  a  chaotic  solution.  The  traditional  techniques  that  were  based  on  the  clas¬ 
sical  signal/  noise  paradigm  must  be  modified  to  process  data  generated  by  processes  in 
which  the  signal  and  noise  cannot  be  so  separated. 

Herein  we  have  shown  that  low-dimensional  deterministic  dynamical  systems  can 
generate  time  series  that  appear  to  have  many  of  the  characteristics  of  random  noise.  We 
can,  however,  distinguish  between  a  chaotic  and  a  completely  random  time  series  by 
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means  of  the  attractor  reconstruction  technique.  The  time  series  data  used  in  the  recon¬ 
struction  of  attractors  can  be  the  vertical  displacement  of  the  sea  surface,  a  component  of 
the  horizontal  fluid  velocity  at  a  fixed  point  on  the  sea  surface  or  in  its  interior,  measure¬ 
ments  at  multiple  points,  and  so  on.  These  data  may  then  be  used  to  reconstruct  the 
dynamic  attractors  determining  the  local  behavior  of  the  sea.  The  underlying  assumption 
is  that  the  nonlinear  hydrodynamic  interactions  are  sufficiently  strong  that  a  single  vari¬ 
able  is  influenced  by  each  of  the  other  variables  to  such  a  degree  that  the  time  lagged 
variables  contain  the  same  information  as  the  original  set  of  variables  characterizing  the 
hydrodynamics.  It  is  worth  mentioning  that  the  more  variables  that  are  measured,  veloci¬ 
ties  at  multiple  space  points,  or  the  vertical  velocity  of  the  surface  as  well  as  its  displace¬ 
ment,  all  give  additional  confidence  in  the  attractor  reconstructed  from  the  time  series. 

Of  course,  the  existence  of  such  attractors,  are  of  interest  to  physical  oceanogra¬ 
phers,  but  what  is  perhaps  of  more  interest  is  how  the  ambient  attractor  could  be  modified 
by  the  disturbances  generated  by  a  surface  current.  Said  differently,  what  we  want  to 
accomplish  is  to  use  the  time  series  X  ( t)  to  predict  the  future  evolution  of  the  wave  field 
beyond  the  available  data  and  to  examine  how  various  currents  disrupt  this  evolution. 

4. 

Barnsley  (1986)  has  been  able  to  use  chaos  to  predict  chaos.  He  developed  an  algorithm' 
to  replicate  a  chaotic  curve  generated  by  a  laboratory  experiment  at  a  finite  number  of 
points.  The  resulting  numerically  generated  curve  “fits”  the  experimental  curve  surpris¬ 
ingly  well  both  inside  and  outside  the  fitted  interval. 

Previous  processing  techniques  assume  the  wave/wave  and  other  hydrodynamic 
interactions  to  be  weak  and  predict  a  weak  modulation  of  the  surface  wave  spectrum  in 
the  Bragg  wave  region.  Such  analyses  have  been  based  on  perturbation  theory,  keeping 
only  the  lowest  order  interactions  between  the  surface  waves  and  surface  currents.  The 
signal  is  then  given  by  the  modulation  in  the  clutter  cross  section  of  the  sea  surface  and  is 
directly  proportional  to  the  surface  current  [Thompson  and  West,  1975].  Improvements 
on  this  approach  involve  integrating  along  the  direction  of  the  current  to  provide  a 
coherent  gain  in  the  signal;  using  multiple  frequency  radars  to  sample  the  modulation  in 
different  spectral  domains  as  well  as  others.  We  do  not  have  room  here  to  critique  the 
relative  merits  of  each  of  these  techniques.  We  can  point  out,  however,  that  none  of 


•  He  regards  the  experimentally  generated  curve  as  the  chaotic  attractor  of  some  descriptive  map.  Then  he  constructs  this 
dissipauvc  map  by  forming  the  union  of  a  sequence  of  dissipative  maps  randomly  selected  from  a  predetermined,  finite  set. 
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them  take  the  basic  nonlinear  hydrodynamic  aspects  of  the  sea  surface  and  interior  into 
account  in  their  processing.  The  attractor  reconstruction  approach  is  the  only  one  that 
intrinsically  depends  on  the  nonlinear  nature  of  the  time  series  for  its  implementation. 

In  the  same  way  the  contamination  of  the  acoustic  signal  by  undersea  phenomena  is 
treated  perturbatively,  relying  on  the  classic  signal/noise  paradigm  to  describe  the 
influence  of  internal  waves  and/or  turbulence  on  the  propagating  pressure  field.  Again 
there  are  established  techniques  similar  to  those  described  above  intended  to  enhance  the 
signal-to-noise  ratio  that  do  not  capitalize  on  the  intrinsic  nonlinearity  of  the  dynamics  of 
the  intervening  medium  and  its  interaction  with  the  acoustic  wave.  It  is  not  yet  clear  how 
these  traditional  techniques  apply  when  the  dynamics  are  chaotic.  We  note  that  there 
already  exists  a  simple  model  of  internal  waves  based  on  the  strange  attractor  concept 
[Abarbanel,  1983]. 

The  relation  between  the  hydrodynamic  patterns  observed  in  nature  and  the  non¬ 
linear  dynamics  of  fluid  particles  has,  over  the  past  decade,  become  increasingly  well 
understood.  Although  fluid  mixing  remains  a  mystery,  some  of  the  concepts  underlying 
the  formation  of  coherent  structures  and  their  eventual  breakup  due  to  chaos  are  becom¬ 
ing  clearer  [Wiggins,  1988;  Ottino,  Lsong,  Rising  and  Swanson,  1988]. 
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Figure  (2.1.2)  The  collection  of  trajectories  initiated  from  a  set  of  initial  conditions  is 
called  a  flow  field.  Here  the  flow  field  in  the  neighborhood  of  a  fixed  point  is  shown.  All 
the  trajectories  asymptotically  converge  on  the  single  point  in  phase  space,  i.e.,  they 
reside  in  the  basin  of  attraction  of  the  focus. 
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Figure  (2.1.3)  The  two  points  a  and  b  in  the  figure  are  possible  initial  conditions  for  the 
system.  When  the  system  can  manifest  limit  cycle  behavior  the  orbits  approach  this 
cycle  asymptotically  and  lose  all  memory  of  their  initial  state. 


Figure  (2.1.4)  The  (x,;c)  phase  space  is  shown  for  a  harmonic  oscillator  with  a  few  typi¬ 
cal  orbits.  Each  ellipse  has  a  constant  energy.  As  the  energy  of  the  oscillator  is 
increased  the  system  jumps  from  an  ellipse  of  smaller  diameter  to  one  of  larger  diameter. 
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Figure  (2.1.6)  The  time  history  of  the  Y  (t)  component  of  the  solution  to  the  Lorenz  sys¬ 
tem  of  equations  C2.1.2)-(2.1.4)  is  shown  for  3xl03  time  units  (from  Lorenz,  1963). 
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Figure  (2.1.10)  The  "funnel”  attractor  solution  to  the  Rossler  Equations  (2.1.12)- 
(2.1.14)  with  parameter  values  a  =0.343, b  ~  1.82  and  c  =9.75.  (From  Rossler.  1979). 
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Figure  (2.1.11)  An  x-y  phase  plane  plot  of  the  solution  to  the  Rossler  Equaitons 
(2.1.1 2)-(2. 1. 14)  with  parameter  values  a  =  0.20  and  b  =  0.20  at  four  different  values 
of  c  indicated  in  the  graphs. 


Figure  (2.2.1)  The  spiral  is  an  arbitrary  orbit  depicting  a  function  y  =f  (X).  The  inter¬ 
section  of  the  spiral  with  the  x-axis  defines  a  set  of  points  x^x^  ■■■,  that  can  be  obtained 
from  a  mapping  determined  by  /  ( x ). 
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vith  a  single  maximum  is  shown.  In  (a),  the  iteration 
:cted.  In  (b),  the  convergence  to  the  station  point  y  * 
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Figure  (2  2  5)  The  bifurcation  of  the  solution  to  the  mapping  1  -  \xx2  as  a  function 
of  ^  -  [L  is  indicated.  The  logarithmic  scale  was  chosen  to  clearly  depict  the  bifurcation 

regions.  (From  Collett  and  Eckmann,  1980). 


Figure  (2.2.7)  An  arbitrary  trajectory  is  shown  and  its  intersection  with  a  plane  parallel 
to  the  xj, jc3  -  plane  at  x2  =  constant  are  recorded.  The  points  A,  B,  C,...  define  a  map  as 
in  Figure  (2.2.1).  This  is  the  Poincard  surface  of  section.  (From  Ott,  1985). 
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Figure  (2.2.9)  (a)  Enlargement  of  the  boxed  region  in  Figure  (2.2.8),  104  iterations;  (b) 
enlargement  of  the  square  in  (a),  106  iterations;  (c)  enlargement  of  the  square  in  (b), 
5  x  106  iterations.  (From  Ott,  1985). 
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Figure  (2.2.10)  Lyapunov  exponents  define  the  average  stretching  or  contraction  of  tra¬ 
jectories  in  characteristic  directions.  Here  we  show  the  effects  of  applying  a  two- 
dimensional  mapping  to  circles  of  initial  conditions.  A  sufficiently  small  circle  of  radius 
e  is  transformed  after  n  iterations  into  an  ellipse  with  major  radius  X."  jE  and  minor  radius 
Xn2e,  where  Xx  and  X2  are  the  Lyapunov  exponents  for  n  -><*>. 


Figure  (2.3.1)  The  time  trace  of  a  random  function  X(r)  versus  time  t  is  shown  in  the 
upper  curve.  The  lower  curve  is  the  same  time  trace  displaced  by  a  time  interval  x.  The 
product  of  these  two  functions  when  averaged  yield  an  estimate  of  the  autocorrelation 
function  (?„( t,T). 


FREQUENCY (u-u0) 

Bgure  (2.3.2)  (a)  The  autocorrelation  function  C«(t)  for  the  typical  time  traces  dep¬ 
icted  m  Figure  (2.3.1)  assuming  the  fluctuations  are  exponentially  con-elated  in  time 
lexpt -x/xc )].  The  constant  xc  is  the  time  required  for  C„ (x)  to  decrease  by  a  factor  1/e, 
is  is  the  decorrelation  time,  (b)  The  power  spectral  density  Sa  (co)  is  graphed  as  a  func¬ 
tion  of  frequency  for  the  exponential  correlation  function  with  a  central  frequency  co0. 
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Figure  (2.3.3)  The  power  spectral  density  for  the  X  (t )  time  series  for  the  “funnel”  dep¬ 
icted  in  Figure  (2.1.10). 
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Figure  (2.3.4)  (a)  The  correlation  integral  for  the  logistic  map  (2.2.11)  at  the  infinite 
bifurcation  point  p  =  |i„  =  3.699  •••  The  starting  point  was  Y0=l/2,  the  number  of 
points  was  N  =3xl04.  (b)  Correlation  integral  for  the  Henon  map  (2.2.29)  and  (2.2.30) 
with  c  - 1.1,13=0.01  and  (V  =  1.5X104.  (c)  Correlation  integrals  for  the  Lorenz  equations 
(2. 1 .9)-(2. 1.1 1)  (dots);  for  the  Rabinovich-Fabricant  equation  (open  circles).  In  both 
cases  A/  =  1.5xl04  and  x  =  0.25.  (From  Grassberger  and  Procaccia,  1985).  (From 
Grassburger  and  Procaccia,  1983a). 
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Figure  (3.1.1)  Time  trajectory  of  three-hourly  surface  pressure  evolving  in  a  three- 
dimensional  phase  space  of  time-lagged  pressure  coordinates  (units:  1000  mb).  Horizon¬ 
tal  axis:  p(t):  vertical  axis:  p(t  +r);  axis  into  plotting  plane:  p(t  +2 r).  From  left  to 
right:  (a)  r  =  3  hours,  (b)  r  =  1  dav,  (c)  r  =  3  days  (from  Fraedrich,  1986). 
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Figure  (3.1.2)  Cumulative  distribution  Junction  of  distances  of  the  15-year  daily  pres¬ 
sure  trajectory  evolving  in  m  -dimensional  phase  spaces  (m  =  1  to  20)  of  time  lagged 
(r  =  3  days)  coordinates  of  the  same  varible.  Top:  observed  time  series;  bottom:  related 
random  series  of  same  length,  mean  and  variance  (from  Fraednch,  1986). 
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Figure  (3.1.3)  Dimensionality  (d)  of  the  weather  attractor  as  a  function  of  the  number  in 
of  phase  space  coordinates,  which  are  determined  by  multiples  ( m  )  the  time  shifted  pres¬ 
sure  record  (using  lags  x  =  1,3  and  6  days).  Randomized  time  series  (index  c)  are  com¬ 
pared  with  observations.  From  left  to  right:  (a)  15  year  record,  (b)  winter  seasons,  (c) 
summer  seasons  (from  Fraedrich,  1986). 
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Figure  (3.1.5)  Summary  of  results  of  scaling  exponent,  d(m)  as  a  function  of  embed¬ 
ding  dimension  m,  obtained  by  two  methods  (method  1,  +;  method  2,  A).  A  random 
number  set  of  the  same  size  as  'he  data  set  was  used  as  a  control  (  □  ).  Note  the  satura¬ 
tion  of  the  exponent  arising  from  the  data,  for  both  methods,  while  there  is  no  saturation 
for  the  random  number  set.  This  limiting  value  is  identified  as  v  (from  Essex,  Lookman 
and  Nerenberg,  1987). 
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Figure  (3.1.7)  (a)  The  data  used  in  this  study  represent  10-second  averages  of  the  verti¬ 
cal  wind  velocity  over  1 1  hours.  The  air  close  to  the  ground  is  heated  and  rises,  creating 
strong  convection.  Positive  values  indicate  updrafts  and  negative  values  indicate  down- 
drafts.  (b)  The  autocorrelation  function  for  the  above  data.  The  inset  graph  is  a 
magnification  of  the  region  close  to  the  origin,  (c)  The  logarithm  of  the  spectral  density 
as  a  function  of  the  frequency  for  the  above  data.  The  spectra  show  various  peaks  on  a 
background  of  a  continuous  frequency  spectrum.  This  suggests  that  a  strange  uiu'ucioi 
may  be  present  (from  Tsonis  and  Eisner,  1988). 


Figure  (3.1.8)  Plot  of  In  Cm(r)  against  In  r  for  embedding  dimensions, 
m  =  4,6,8,10, 12.  Note  the  convergence  of  slopes  as  m  increases  (from  Tsonis  and 
Eisner,  1988). 


Figure  (3.1.9)  Scaling  exponent,  d(m),  as  a  function  of  the  embedding  dimension,  m. 
Crosses  correspond  to  the  wind  velocity  data  and  squares  to  a  random  sample  of  the  same 
size  as  the  wind  data.  Note  the  saturation  of  the  scaling  exponent  observed  for  the  wind 
data,  although  there  is  no  saturation  for  the  random  set.  From  this  figure  we  estimated 
v  =  (from  Tsonis  and  Eisner,  1988). 


Figure  (3.2.4)  Four  magnifications  of  the  surface  for  D  =  2.05,  showing  self-similarity 

( M  =8 ,b  =  1.2).  The  upper  right  surface  is  the  fivefold  magnification  of  a  section  of 

the  upper  left  surface.  Similarly,  the  lower  left  surface  is  a  fivefold  magnification  of  a  •« 

piece  of  the  upper  right  surface  and  the  lower  right  surface  is  a  fivefold  magnification  of 

the  lower  left  surface.  The  vertical  extent  is  magnified  by  (from  Ausloos  and 

Berman,  1985). 
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Figure  (4.1.2)  The  correlation  dimension  d{m)  is  plotted  versus  the  embedding  dimen¬ 
sion  m  for  the  data  in  Figure  (4.1.1).  The  solid  line  depicts  the  correlation  d  (m )  =  m 
for  a  calculated  random  surface.  The  correlation  dimension  are  as  shown  for  a  fully 
developed  one-dimensional  sea  surface  using  the  “time"  shifts  x  =  8(*)  and  t  =  15(x). 
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Figure  (4.1.3)  Standard  correlation  integrals  for  autocorrelated  random  data,  (a)  uB- 
log  plot  of  the  standard  (lV  =  1)  correlation  integral  over  a  range  of  embedding  dimen¬ 
sion  m  for  stochastic  data  with  N  =  10000  points,  standard  deviation  0  =  20,  and  auto¬ 
correlation  a  =  0.9. 
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Figure  (4.1.5)  The  correlation  dimension  d(m)  is  plotted  versus  the  embedding  dimen¬ 
sion  m  for  the  data  in  Figure  (4.1.4).  The  time  series  clearly  indicates  a  saturation  in  the 
neighborhood  of  v  =  6. 


